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Abstract 

Robust  parameter  design  (RPD)  is  used  to  identify  a  system’s  control  settings  that 
offer  a  compromise  between  obtaining  desired  mean  responses  and  minimizing  the 
variability  about  those  responses.  Two  popular  combined-array  strategies — the  response 
surface  model  (RSM)  approach  and  the  emulator  approach — are  limited  when  applied  to 
simulations.  In  the  former  case,  the  mean  and  variance  models  can  be  inadequate  due  to 
a  high  level  of  non-linearity  within  many  simulations.  In  the  latter  case,  precise  mean 
and  variance  approximations  are  developed  at  the  expense  of  extensive  Monte  Carlo 
sampling. 

This  research  combines  the  RSM  approach’s  efficiency  with  the  emulator 
approach’s  accuracy.  Non-linear  metamodeling  extensions,  namely  through  Kriging  and 
radial  basis  function  neural  networks,  are  made  to  the  RSM  approach.  The  mean  and 
variance  of  second-order  Taylor  series  approximations  of  these  metamodels  are  generated 
via  the  Multivariate  Delta  Method  and  subsequent  optimization  problems  employing 
these  approximations  are  solved.  Results  show  that  improved  prediction  models  can  be 
attained  through  the  proposed  approach  at  a  reduced  computational  cost.  Additionally,  a 
multi-response  RPD  problem  solving  technique  based  on  desirability  functions  is 
presented  to  produce  a  solution  that  is  mutually  robust  across  all  responses.  Lastly, 
quality  measures  are  developed  to  provide  a  holistic  assessment  of  several  competing 
RPD  strategies. 
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NON-LINEAR  METAMODELING  EXTENSIONS  TO  THE 


ROBUST  PARAMETER  DESIGN  OF  COMPUTER  SIMULATIONS 


I.  Introduction 


1.1  Background 

Computer  simulations  are  mathematical  models  of  complex  real-world  systems 
for  which  it  would  be  too  expensive,  too  time  consuming,  too  dangerous,  or  even 
impossible  to  examine  with  physical  experimentation.  They  have  been  applied  across 
such  diverse  disciplines  as  manufacturing,  military  applications,  logistics  and  supply 
chain  management,  and  health  care.  Simulations,  which  are  either  deterministic  or 
stochastic  in  nature,  can  be  used  to  gain  insight  into  a  system,  examine  “what-if  ’ 
scenarios,  or  perform  system  optimization  [1].  Detenninistic  simulations  yield  the  same 
outputs  when  multiple  replications  are  perfonned  at  a  single  design  setting.  Stochastic 
simulations,  on  the  other  hand,  generate  different  outputs  when  replications  are 
perfonned  at  a  single  design  setting  due  to  the  presence  of  random  variables  within  the 
model  [2], 

This  research  focuses  on  stochastic  simulations  in  which  some  parameters,  called 
control  factors,  can  be  set  to  specific  values  whereas  other  parameters,  called  noise 
factors,  are  modeled  by  random  variables  with  specific  probability  distributions.  The 
noise  factors  introduce  undesirable  variation  in  the  system’s  outputs.  When  noise  factors 
are  present,  robust  parameter  design  (RPD)  principles  can  be  applied  as  a  cost-effective 
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strategy  for  identifying  the  ideal  control  factor  settings  that  result  in  outputs  that  are 
robust,  or  insensitive,  to  the  random  fluctuation  of  the  noise  factors.  Simulations  are 
well-suited  for  RPD  studies  since  the  parameters,  whether  they  are  control  factors  or 
noise  factors,  can  be  manipulated  and  set  to  specific  values.  In  fact,  Kleijnen  et  al.  [3] 
state  that  finding  robust  policies  or  decisions  is  one  of  the  basic  goals  of  simulations. 
RPD  is  a  method  for  detennining  the  control  factor  settings  that  reach  a  compromise 
between  obtaining  a  desired  mean  response  while  minimizing  the  variability  about  that 
response  [4]. 

In  the  single-response  RPD  problem,  the  objective  is  to  identify  the  control  factor 
setting  that  yields  a  desired  mean  response  with  minimum  variance  [4],  The  literature 
offers  several  strategies  that  have  stemmed  from  Taguchi’s  original  RPD  principles. 
However,  this  research  is  motivated  by  those  strategies  that  utilize  Welch  et  al.’s  [5] 
combined-array  design  of  experiment  (DOE).  Specifically,  the  methods  of  interest  are 
the  combined-array  response  surface  model  (RSM)  approach  [6-9]  and  the  stochastic 
emulator  approach  [10-13].  Though  each  strategy  has  its  own  merits,  they  can  be 
inappropriate  or  cumbersome  when  faced  with  the  highly  non-linear  nature  of  typical 
simulations. 

In  the  multi-response  RPD  problem,  the  objective  is  to  find  the  optimal  control 
parameter  levels  that  return  average  responses  close  to  their  target  values  while 
minimizing  the  variance  of  each  response.  The  literature  proposes  several  methods  for 
optimizing  multi -response  problems.  These  methods  involve  the  use  of  desirability 
functions  [14-18],  loss  functions  [19-23],  principal  component  analysis  (PCA)  [24-27], 
distance  metrics  [28,  29],  and  mean  square  error  (MSE)  criterion  [30-32].  A  majority  of 
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these  techniques  transform  the  set  of  quality  characteristics  into  new  response  variables  in 
order  to  reduce  the  dimension  of  the  optimization  problem.  These  existing  methods  seek 
to  find  an  optimal  balance  of  means  and  variances  across  the  set  of  responses.  However, 
in  some  instances,  the  mean  or  variance  of  one  response  may  influence  the  solution  in 
such  a  way  that  the  means  and  variances  of  the  remaining  responses  are  insignificant  to 
the  overall  RPD  problem.  In  this  case,  it  can  be  difficult  to  attain  a  solution  that  is 
balanced  across  the  set  of  responses. 

Numerous  organizations  utilize  RPD  principles  as  an  economical  strategy  for 
developing  a  product  or  process  that  is  insensitive  to  a  variety  of  operating  conditions. 
Advantages  of  employing  these  principles  are  that  the  product  or  process  will  be  on 
target  while  exhibiting  less  variability.  This  subsequently  increases  the  end  user’s  appeal 
for  the  product  or  process  since  it  won’t  be  as  susceptible  to  deterioration  and  can  be  used 
in  diverse  situations.  This  research  effort  aims  to  advance  the  current  RPD  problem 
solving  approaches. 

1.2  Research  Objectives 

This  dissertation  strives  to  meet  three  key  objectives.  The  first  objective  is  to 
broaden  the  combined-array  RSM  approach  that  relies  exclusively  on  low-order 
polynomial  models.  Since  more  accurate  predictive  response  surface  models  result  in 
better  RPD  solutions  [33],  a  methodology  will  be  developed  that  utilizes  non-linear 
modeling  efforts,  such  as  Kriging  and  radial  basis  function  neural  networks  (RBFNNs), 
in  place  of  polynomial  models.  The  second  objective  is  to  develop  an  approach  for  multi¬ 
response  RPD  problems  that  provides  a  collaborative  solution  that  is  balanced  across  the 
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means  and  variances  of  each  response.  Finally,  the  third  objective  is  to  generate  a 
framework  for  comparing  different  RPD  problem  solving  strategies  via  quality  measures. 
Such  measures  can  increase  the  understanding  of  each  approach  and  allow  the  analyst  to 
make  a  more  knowledgeable  evaluation  of  the  competing  procedures. 

1.3  Chapter  Overview 

This  dissertation  is  organized  in  the  following  manner.  Chapter  II  establishes  the 
foundation  for  the  techniques  utilized  in  this  research  by  reviewing  the  pertinent  literature 
in  the  areas  of  RPD  and  response  modeling.  Chapter  III  details  the  extension  of  the 
combined-array  RSM  approach  to  include  the  application  of  Kriging  and  RBFNN 
metamodels.  Chapter  IV  proposes  a  multi-response  RPD  methodology  based  on 
desirability  functions  that  generates  solutions  that  are  well-balanced  across  the  means  and 
variances  of  each  response.  Chapter  V  describes  a  method  for  comparing  RPD  problem 
solving  strategies  via  quality  measures.  Chapter  VI  summarizes  the  research 
contributions  and  identifies  opportunities  for  future  research.  Finally,  supplemental 
mathematical  derivations  and  figures  are  presented  in  the  Appendices. 
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II.  Pertinent  Literature 


2.1  Introduction 

This  chapter  summarizes  pertinent  literature  in  the  areas  of  RPD  and  response 
modeling.  It  is  intended  to  establish  the  foundation  for  the  techniques  that  are  utilized 
and  expanded  upon  in  this  dissertation. 

2.2  Robust  Parameter  Design  Approaches 

This  research  focuses  on  the  robust  design  of  stochastic  simulations  in  which 
some  parameters,  called  control  factors,  can  be  set  to  specific  values  whereas  other 
parameters,  called  noise  factors,  are  modeled  by  random  variables  with  specific 
probability  distributions.  The  noise  factors  introduce  undesirable  variation  in  the 
system’s  output.  When  noise  factors  are  present,  RPD  principles  can  be  applied  as  a 
cost-effective  strategy  for  identifying  the  ideal  control  factor  settings  that  result  in  outputs 
that  are  robust,  or  insensitive,  to  the  random  fluctuation  of  the  noise  factors. 

Simulations  are  well-suited  for  RPD  studies  since  the  parameters,  whether  they 
are  control  factors  or  noise  factors,  can  be  manipulated  and  set  to  specific  values.  RPD  is 
a  method  for  determining  a  system’s  ideal  control  factor  setting  that  reaches  a 
compromise  between  obtaining  a  desired  target  on  the  mean  response  while  minimizing 
the  variability  of  the  system  around  that  mean  response.  Three  popular  strategies  for 
solving  RPD  problems  are  Taguchi’s  method  [4],  the  response  surface  model  (RSM) 
approach  [6],  and  the  stochastic  emulator  approach  [11].  However,  before  discussing 
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these  current  strategies,  this  section  will  first  illustrate  the  sources  that  cause  variation  in 
a  system’s  output. 

2.2.1  Sources  of  Output  Variation 

As  seen  in  Figure  la,  a  system  can  be  modeled  as  a  set  of  resources  that 
transforms  inputs  into  measureable  outputs.  This  figure  expands  Figure  1-1  in 
Montgomery  [9]  to  simulations.  Noise  causes  variability  in  the  outputs  through  any 
combination  of  four  separate  sources:  uncontrollable  factors,  uncertainty  in  the  control 
factors,  uncertainty  in  the  inputs,  and  internal  system  variability.  Uncontrollable  factors 
are  system  parameters  that  are  difficult,  costly,  or  impossible  to  control  under  nonnal 
operation  of  the  system.  In  industrial  scenarios,  these  factors  are  usually  assumed 
controllable  for  the  purposes  of  experimentation.  Uncertainty  in  a  control  factor  reflects 
variation  about  its  nominal  setting.  Uncertainty  in  an  input  involves  variability  of  the 
system’s  inputs.  Finally,  internal  variability  relates  to  intrinsic  system  variation.  For 
example,  in  a  discrete  event  simulation,  intrinsic  system  variation  occurs  due  to  factors 
such  as  the  initial  state  of  the  system,  the  warm-up  period,  the  tennination  conditions,  or 
the  random  number  stream. 

In  this  research,  the  term  noise  factor  represents  uncontrollable  factors, 
uncertainty  in  the  control  factors,  and  uncertainty  in  the  inputs.  The  set  of  noise  factors  is 
denoted  by  the  vector  z  whereas  the  set  of  control  factors  is  specified  by  the  vector  x. 
Typically,  it  is  assumed  that  the  noise  factors  are  mutually  independent  and  that  each 
noise  factor  z;  is  a  normally  distributed  random  variable  with  known  mean  jdi  and  variance 
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cr2  .  Therefore, E[z]  =  pz  =  [//j  •  •  •  jun]  and  Var  [z]  =  Ez  =  diag(crf , of  ) .  Figure 
lb  illustrates  the  sources  of  variation  for  a  scenario  used  in  Chapter  III. 
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Figure  1,  Illustrations  of  Sources  of  Output  Variation  for  (a)  a  General  System  and 

(b)  a  Circuit  Simulation 


2.2.2  Taguchi’s  Method 

Taguchi’s  approach  to  solving  RPD  problems  crosses  an  inner  orthogonal  array  of 
control  factors  with  an  outer  orthogonal  array  of  noise  factors.  That  is,  each  combination 
of  control  factor  settings  within  the  inner  array  is  perfonned  over  each  combination  of 
noise  factor  settings  in  the  outer  array.  Taguchi  then  summarizes  the  observations  of 
each  inner  array  trial  across  the  outer  array  settings  via  a  statistic  known  as  a  signal-to- 
noise  ratio  (SNR).  The  SNR  accounts  for  both  the  process  mean  and  variance.  It  is  then 
utilized  as  the  response  variable  for  statistical  analysis  purposes  [4], 

The  goal  of  the  overall  experiment  detennines  which  of  four  SNR  formulations  to 
use.  If  the  goal  is  to  determine  the  values  of  the  control  factors  that  result  in  a  minimum 
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response,  then  the  smaller-the-better  SNR  should  be  utilized.  If  a  maximum  response  is 
desired,  then  the  larger-the-better  SNR  is  used.  Finally,  in  what  is  known  as  the  target- 
is-best  case,  the  goal  may  be  to  detennine  the  control  settings  that  yield  a  response  near 
some  desired  target  value.  In  this  case,  Taguchi  first  recommends  using  bias-eliminating 
tuning  factors  that  result  in  an  expected  process  response  equal  to  the  target.  Then  the 
form  of  the  SNR  is  based  on  whether  or  not  the  response  mean  and  variance  are 
independent  [9,  34].  These  four  SNRs  are  shown  in  Table  1  where  yj  .  is  the  response  for 

inner  array  setting  i  and  outer  array  settingy,  yt  is  the  sample  mean  of  inner  array  setting 

i,  and.sf  is  the  sample  variance  of  inner  array  setting  i. 


Table  1.  Taguchi’s  Signal-to-Noise  Ratios  for  Three  Experimental  Scenarios 


Smaller-the-Better 

Larger-the-Better 

Target-is-Best 

Independent  mean  &  variance 

(  «  V2  > 

SNR,  =  10 -log  V  l’J 

O-i  n ) 

f  "  l/v2 .  ^ 

SNR,  =  10 -log  V  /  J 

0-  n  j 

SNR, .  =10-log(j,2) 

Correlated  mean  &  variance 

SNR,  =  1 0  ■  log  \ 

\  s<  J 

Regardless  of  which  SNR  is  utilized,  modeling  analysis  of  the  system  attempts  to 
maximize  the  SNR  while  driving  the  mean  response  towards  some  preferred  target  value. 
Taguchi  methods  typically  develop  only  main-effects  models  and  are  not  concerned  with 
control  factor  interactions.  Analysis  of  Taguchi’s  approach  is  a  two-stage  method.  First, 
the  experimenter  should  choose  the  levels  of  those  significant  control  factors  that 
maximize  the  SNR.  Second,  the  experimenter  should  choose  the  levels  of  the  remaining 
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control  factors  that  result  in  a  mean  response  near  the  desired  target  value  [34],  Standard 
ANOVA  techniques  can  also  help  distinguish  which  control  factors  affect  the  average 
response  and  SNR.  Control  settings  that  are  robust  and  insensitive  to  the  variance  caused 
by  the  noise  factors  are  then  detennined  [9]. 

If  it  can  be  assumed  that  there  are  no  significant  interactions  among  the  control 
variables,  then  Taguchi’s  methodology  is  highly  appropriate  for  identifying  robust 
control  parameter  settings.  However,  Taguchi’s  critics  claim  that,  more  often  than  not, 
this  assumption  does  not  hold  and  that  a  main  effects-only  study  may  yield  ambiguous 
results  [7,  35].  Detractors  also  note  that  the  use  of  Taguchi’s  SNRs,  though  they  are 
concerned  with  the  process  mean  and  variance,  don’t  allow  for  a  complete  understanding 
of  which  control  factors  affect  the  mean  and  which  affect  the  variance  [9,  34,  36].  Nair 
and  Shoemaker  [37]  further  argue  that  critical  system  information  is  lost  by  compressing 
the  experimental  responses  into  SNRs.  Several  authors  claim  that  it  is  better  to  examine 
the  mean  and  variance  separately  instead  of  combining  them  into  a  single  SNR  [38-40]. 
Even  though  the  purpose  behind  SNRs  is  to  uncouple  the  location  and  dispersion  effects, 
Montgomery  [9]  contends  that  there  is  no  assurance  that  this  will  occur.  As  an  example, 
he  shows  how  the  use  of  the  smaller-the-better  SNR  actually  confounds  location  and 
dispersion  effects.  Pignatiello  [41]  states  that,  despite  the  criticism  of  Taguchi’s  tactics, 
his  conceptual  framework  for  planning  a  product  or  process  design  experiment  is 
fundamentally  sound. 


9 


2.2.3  Response  Surface  Model  Approaches 

The  response  surface  model  (RSM)  approach  solves  an  optimization  problem 
involving  models  of  the  system’s  mean  and  variance.  Several  optimization  schemes  are 
discussed  in  Section  2.4.  Two  major  RSM  applications  are  found  in  the  literature.  One 
approach  utilizes  replications  of  an  experimental  design  whereas  the  second  approach 
uses  a  combined-array  experimental  design. 

2.2. 3.1  RSM  Approach  Using  Replicated  Experimental  Designs 

The  first  RSM  application  uses  replications  of  a  design  of  experiment  (DOE) 
consisting  of  x  only.  This  approach  is  used  when  the  system’s  output  variability  is  due  to 
either  random  sampling  from  z  or  intrinsic  system  variation.  Separate  low-order 
polynomials  [6],  higher-order  polynomials  [33,  42],  Kriging  models  [11,  13,  43],  radial 
basis  function  approximations  [43],  or  radial  basis  function  neural  networks  [43]  are  fit  to 
the  sample  means  and  variances  of  the  design  points  to  generate  mean  and  variance 
models.  This  RSM  application  has  been  applied  to  simulations  of  a  piston  [11]  and  an 
optical  profilometer  [42].  Wild  and  Pignatiello  [44]  used  a  partitioning  strategy  to  assess 
the  system’s  mean  and  variance  through  crossed-array  designs.  They  used  a  discrete 
event  simulation  for  the  robust  design  of  a  manufacturing  facility.  Several  authors  have 
also  applied  this  RSM  application  to  Box  and  Draper’s  [45]  well-examined  printing 
process  study  [33,  43,  46-48]. 

2. 2. 3. 2  RSM  Approach  Using  Combined-Array  Experimental  Designs 

The  second  RSM  application  uses  the  combined-array  DOE  proposed  by  Welch  et 
al.  [5].  The  combined-array  is  comprised  of  both  x  and  z.  A  low-order  polynomial 
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response  model  j)(x,  z)  is  built  to  accommodate  the  linear  control  factor  and  noise  factor 
effects,  the  pure  quadratic  control  factor  effects,  the  control  factor  interaction  effects,  and 
the  control  factor  by  noise  factor  interaction  effects  [7].  Mean  and  variance  models  are 
then  obtained  by  taking  the  expectation  and  variance  of  j)(x,  z) .  Recent  extensions  have 
been  made  to  include  pure  quadratic  noise  factor  effects  and  noise  factor  interaction 
effects  [49]  as  well  as  three-factor  control  by  control  by  noise  factor  interactions  and 
control  by  noise  by  noise  factor  interactions  [50].  The  combined-array  RSM  approach 
has  been  applied  to  a  piston  simulation  [11]  and  an  economic  order  quantity  inventory 
model  [12,  13].  A  textbook  example  of  this  approach  for  a  physical  experiment  can  be 
found  for  a  chemical  production  process  in  Montgomery  [9]. 

2.2.4  Stochastic  Emulator  Approaches 

The  literature  discusses  two  main  applications  of  the  stochastic  emulator  approach 
for  solving  RPD  problems  for  stochastic  simulations.  Regardless  of  the  application,  the 
first  step  is  to  generate  a  model,  or  emulator,  of  the  system.  The  emulator  is  then  used  to 
generate  large  samples  of  data  from  the  joint  distribution  of  z  at  the  design  points  in  x. 
Mean  and  variance  models  are  fit  from  the  samples  and  a  subsequent  optimization 
problem  consisting  of  the  mean  and  variance  models  is  solved. 

2.2. 4.1  Emulator  Approach  with  Control  Factor  Uncertainty 

In  the  first  emulator  application,  uncertainty  in  the  control  factors  causes  variation 
in  the  system’s  output.  In  this  case,  control  factor  setting  xt  can  be  decomposed  into  the 
sum  of  the  nominal  control  factor  setting  x.  and  the  normally  distributed  noise  factor  z. . 
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That  is,  X;  =  xt  +  z;  where  A  [x,]  =  xt and  Var  [x; ]  =  a]  ;  hence A[x]  =  x  and Var\x\  =  Zz . 
Initially,  the  noise  factors  are  ignored  by  setting  z  =  0  and  the  emulator  j>(x  =  x)  is  built 
from  a  DOE  in  the  control  factors  x.  A  second  DOE  in  x  is  then  constructed  and  the 
emulator  is  used  to  evaluate  the  noise  factors’  effect  on  the  output.  For  each  design  point 
xd  ,  a  large  sample  of  responses  is  generated  from  the  nonnal  distribution  with  mean  xd 

and  covariance  Zz  as  evaluated  using  the  emulator.  Separate  models  are  fit  to  the  sample 

means  and  variances  of  these  responses.  A  robust  control  setting  is  detennined  by 
optimizing  a  problem  involving  the  emulator’s  mean  and  variance  models.  Bates  et  al. 
[11]  used  this  emulator  strategy  in  an  RPD  study  on  a  piston  simulation. 

2.2. 4.2  Emulator  Approach  with  Uncontrollable  Factors  or  Input  Uncertainty 
In  the  second  emulator  application,  variation  in  the  system’s  output  is  due  to 
uncertainty  in  the  inputs  and/or  the  existence  of  uncontrollable  factors.  The  emulator  is 
built  from  a  combined-array  DOE  including  x  and  z.  The  high  and  low  experimental 
levels  for  each  z;  are  set  at  /ui  ±3cr.  [13],  Once  the  emulator  j)(x,  z)  is  built,  a  second  DOE 

in  x  is  constructed.  For  each  design  point  xd  ,  a  large  sample  is  generated  from  the  joint 

distribution  of  z  and  is  evaluated  using  the  emulator.  As  in  the  first  emulator  approach, 
mean  and  variance  models  are  fit  to  the  sample  means  and  variances  calculated  at  each 
design  point.  These  models  are  then  used  in  an  optimization  formulation  to  find  a  robust 
control  setting.  This  RPD  emulator  strategy  has  been  performed  on  simulations  for  a  2- 
bar  truss  design  problem  [10]  and  an  economic  order  quantity  inventory  model  [13]. 
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2.3  Robust  Parameter  Design  Modeling  Strategies 

The  RPD  approach  developed  in  this  research  will  be  compared  to  the  combined- 
array  RSM  approach  and  the  combined-array  stochastic  emulator  approach.  Thus,  this 
section  will  discuss  the  modeling  efforts  that  have  been  utilized  for  those  two  approaches. 

2.3.1  Low-Order  Polynomial  Models 

There  are  four  modeling  strategies  that  have  been  applied  to  the  combined-array 
RSM  approach.  Each  builds  low-order  polynomial  response  models  comprised  ofrx 
control  factors  x  and/;  noise  factors  z.  Each  approach  discussed  in  this  subsection, 
expands  upon  the  model  that  precedes  it.  They  each  assume  that  the  model’s  error  e  is 
normally  distributed  with  a  mean  of  zero  and  a  constant  variance  a2  ;  hence  any  non¬ 
constant  variance  stemming  from  the  process  is  attributed  to  the  inability  to  control  the 
noise  factors  [8],  They  also  assume  that  each  noise  factor  z;  is  a  normally  distributed 

random  variable  with  a  known  mean  //;  and  variance  of  .  For  experimental  purposes,  the 
natural  level  of  each  noise  factor  is  centered  at  its  mean  and  its  ±1  levels  are  set  at  //,  ±  cr 
.  Furthermore,  the  noise  variables  are  assumed  to  be  uncorrelated.  Therefore,  the  mean 
and  variance  of  z  are  defined  asE’fz]  =  pz  and  Far  [z]  =  Zz  =  diag^af  j  , 

respectively.  Thus,  if  the  factors  have  been  transfonned  to  the  coded  variable  space, 

E\x\  =  0,.  andFar[z]  =  Ir  [9]. 
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2. 3.1.1  The  Standard  RPD  Response  Model 

The  standard  (Std)  fonn  of  the  RPD  response  model  in  Equation  (1)  considers  the 
noise  factor  interactions  and  the  pure  quadratic  noise  factor  effects  negligible  [8]. 

ystd  (x,z)  =  /?0  +  x'p  +  x'Bx  +  (y'  +  x'A)  z  +  z  (1) 

In  Equation  (1),  J30  represents  the  intercept,  x  is  the  rx  x  1  vector  of  control  factors,  z  is  the 
rz  xl  vector  of  noise  factors,  p  is  the  rx  x  1  vector  of  linear  control  factor  effects,  B  is  the 
rx  x  rx  matrix  where  the  pure  quadratic  control  factor  effects  are  on  the  diagonal  and  one- 
half  of  the  control  factor  interaction  effects  are  on  the  off-diagonal,  y  is  the  rz  x  1  vector  of 
linear  noise  factor  effects,  and  A  is  the  rx  x  r2  matrix  of  control  factor  by  noise  factor 
interaction  effects. 

Given  the  distributional  assumptions  of  the  noise  variables  and  the  model’s  error, 
the  mean  and  variance  of  the  standard  response  model  in  Equation  (1)  can  be  written 
respectively  as 

E[yStd  (x,  z)]  =  A,  +  x'p  +  x'Bx  (2) 

and 

Var[ys,d  (x,z)]=  (y'  +  x'A) Zz  (y'  +  x'A)'  +  cr  (3) 

where  cr  is  estimated  using  the  MSE  of  the  fitted  response  surface  model.  Note  that 
Equations  (2)  and  (3)  are  completely  in  terms  of  the  control  factors  x.  These  two 
response  models  can  now  be  used  to  estimate  the  mean  and  variance  of  the  process  at  any 
point  within  the  control  design  space. 
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2.3. 1.2  The  Noise-by-Noise  Response  Model 

Mindrup  et  al.  [49]  found  that  Equation  (1)  was  inadequate  in  certain  imaging 
applications.  As  a  result,  they  removed  the  assumption  that  the  noise  factor  interactions 
and  the  pure  quadratic  noise  factor  effects  are  insignificant  and  developed  the  noise-by- 
noise  (AW)  response  model,  mean  model,  and  variance  model  in  Equations  (4)-(6). 

yNN(x,z)  =  p0  +  x'p  +  x'Bx +  (y'  + x'A)  z  +  z'0>z  +  £  (4) 

E\_ym  (x,  z  )]=/?„+  x'p  +  x'Bx  +  tr  (<PSz )  (5) 

Var[yNN  (x,  z)]  =  ( y'  +  x'A)  Zz  (y'  +  x'A)'  +  2 tr  (<t>Zz<DSz )  +  cr2  (6) 

In  Equations  (4)-(6),  tr  represents  the  trace  of  a  square  matrix  and  <t>  is  therz  xrz  matrix 

with  pure  quadratic  noise  factor  effects  on  the  diagonal  and  one-half  of  the  noise  factor 
interaction  effects  on  the  off-diagonal. 

2. 3. 1.3  The  Control-by-Noise-by-Noise  Response  Model 

Williams  et  al.  [50]  further  extended  Mindrup  et  al.’s  work  in  two  successive 
expansions.  First,  they  appended  the  three-factor  control-by-noise-by-noise  ( CNN) 
interactions  to  generate  the  CNN  RPD  response  model,  mean  model,  and  variance  model 
in  Equations  (7)-(9). 

yCNN  (x,z)  =  /?„  +  x'p  +  x'Bx +  (y'  +  x'A) z  +  z'(®  +  ,Fx)z  +  £  (7) 

E [. yCm  (x> z)]  =  Ao  +  X'P  +  x'Bx  +  tr  ((O  +  Tx )  Iz )  (8) 

Var  [vcw  (x,  z)]  =  (y'  +  x'A) Zz  (y'  +  x'A)'  +  2tr  ((<D  +  XVX )  Zz )"  +  oy2  (9) 
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The  additional  term  in  Equations  (7)-(9)  is  ¥x  =  ^  where  T  is  the  r,  x  rz  matrix  of 

i= 1 

control-by-noise-by-noise  interaction  effects  corresponding  to  control  factor  xt . 

2. 3.1.4  The  Control-by-Control-by-Noise  Response  Model 

Williams  et  al.  [50]  then  further  extended  the  CNN  RPD  models  to  include  the 
three-factor  control-by-control-by-noise  ( CCN)  interaction  terms.  The  CCN  RPD 
response  model,  mean  model,  and  variance  model  are  shown  in  Equations  (10)— (12). 

yCCN  (x,z)  =  J30  +x'p  +  x'Bx  +  (y'  +  x'A  +  x)z  +  z'(0-l-vPx)z  +  ^  (10) 

E  [ yCCN  (x> z)]  =  A,  +  X'P  +  x'Bx  +  flr  ((®  +  Tx )  Sz )  (11) 

^[Tccv(x’z)]  =  (Y'  +  x'A  +  x)Iz(y'  +  x'A  +  x)'+2tr((0  +  Tx)Zz)2+o-£2  (12) 

The  additional  term  in  Equations  ( 1 0)— (12)  is  the  vectorx  =  j^x'QjX,  x'£l2x,  ...  ,  x'Qr  xj 
where  fl .  is  the  r  x  r  matrix  of  control-by-control-by-noise  interaction  effects 
corresponding  to  noise  factor  z . . 

2.3.2  Kriging 

Kriging  has  been  used  to  model  the  simulation’s  response,  as  well  as  its  mean  and 
variance,  within  the  combined-array  stochastic  emulator  strategy  in  Section  2.2.4. 

Kriging  is  a  nonparametric,  global,  exact  interpolation  model.  The  term  nonparametric 
implies  that  training  points  are  used  to  both  estimate  the  unknown  model  parameters  and 
predict  the  response  of  new  observations  [43].  Global  means  that  a  Kriging  model 
provides  predictions  across  the  entire  experimental  area  and  exact  indicates  that  Kriging 
models  predict  the  precise  response  of  previously  observed  input  combinations  [51]. 
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These  properties  are  the  reason  as  to  why  Kriging  metamodels  have  become  popular 
approaches  for  estimating  the  output  of  computer  simulation  models  [2,  11,  43,  5 1-58], 
The  origins  of  Kriging  are  in  geostatistics,  or  spatial  statistics  [59]. 

The  Kriging  model  for  the  input  vector  v  assumes  the  form 

v(v)  =  /(v)  +  £(v).  (13) 

In  Equation  ( 1 3),  /(v)  models  the  trend  in  the  data  and  provides  a  global  approximation  of 
the  design  space  [56].  Ordinary  Kriging  assumes  that  /(v)  =  f]  is  the  constant  mean  of 
the  data  in  the  experimental  region  of  interest  [2]  whereas  Universal  Kriging  uses  a  low- 
order  polynomial  to  define  f{\)  [43].  Also  in  Equation  (13),  S(\)  creates  localized 
deviations  [56]  and  is  additive  noise  fonned  by  a  stationary  covariance  process  with  zero 
mean,  variance  cr2 ,  and  covariance  <t2R  in  which  the  i/h  element  of  the  NxN  spatial 
correlation  matrix  R  is  defined  by 


1 

exp 

vf -v^r 

k= 1 

for  i  =  j 
for  i  ^  j 


(14) 


where  v<k‘  >  is  the  klh  feature  of  the  /th  training  vector,  0k  >  0  is  the  weight  factor  for  the  klh 

input  vector  feature,  p  is  a  parameter  that  defines  the  correlation  between  two  training 
points,  N  is  the  number  of  training  vectors,  and  K  is  the  number  of  features  in  each 
training  vector.  This  research  uses  Ordinary  Kriging  since  that  has  been  found  to  be 
sufficient  for  simulation  models  [43,  54,  60].  The  Gaussian  correlation  function,  where 
p  =  2  ,  is  also  utilized  due  to  its  widespread  appeal  [56]. 
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Estimates  for  the  parameters  /?,  cr2 ,  and  0k  are  detennined  through  the  maximized 


likelihood  estimation  (MLE)  approach.  Interested  readers  can  reference  [43]  for  the 
MLE  derivations.  As  a  result,  given  a  set  of  training  vectors,  the  Kriging  model  response 
for  a  new  observation  v  is 

v(v)  =  /?  +  r'(v)R~1(y  -1/?)  (15) 


-  1R  y 

where  ft  =  — — —is  the  mean  parameter,  y  is  the  N  x  1  vector  of  training  set  responses,  1  is 

1  R  1 


a  N  x  1  vector  of  ones,  and  r(v)  is  the  Nxl  correlation  vector  whose  nth  element  is  the 


correlation  between  v  and  the  nlh  training  vector  v(,,)  defined  as 


,0) , 


r„(v)  =  exp 


k= 1 


(16) 


2.3.3  Radial  Basis  Function  Neural  Networks 

Artificial  neural  networks  (ANNs),  and  in  particular  radial  basis  function  neural 
networks  (RBFNNs),  have  also  been  used  to  provide  a  metamodel  of  a  simulation’s 
response,  mean,  and  variance  within  the  stochastic  emulator  approach.  ANNs  are 
mathematical  models  that  update  their  parameters  iteratively  to  leam  the  relationship 
between  a  set  of  inputs  and  a  set  of  outputs.  The  RBFNN  is  a  special  class  of  ANN  [61- 
64], 

In  a  typical  RBFNN  framework,  as  illustrated  in  Figure  2,  the  number  of  input 
layer  neurons  is  equal  to  the  number  of  input  features  and  the  number  of  hidden  layer 
neurons  is  equal  to  the  number  of  training  vectors.  There  is  also  an  output  layer  neuron 
for  each  system  output.  In  this  case,  there  are  N  training  vectors  with  K  features  and  a 
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single  output  y.  The  activation  function  of  the  nth  neuron,  hn  (•) ,  is  a  radial  basis  function 
designed  to  fire  high  when  a  training  vector  is  very  close  to  the  center  vector  p(I  and  give 
a  diminishing  response  as  the  training  vector  moves  away  from  the  center  vector’s 
receptive  field  defined  by  the  spread  parameter  on .  This  research  employs  the  Gaussian 
response  function 


/?„(v)  =  exp 


'-(v-tOV-i,.)' 

2cr2 


(17) 


The  output  layer  simply  computes  a  linear  weighted  sum  of  the  hidden  layer  neuron 
activations. 


Hi  <*\ 


Flow  Through  Network 


Figure  2.  Structure  of  a  Radial  Basis  Function  Neural  Network 

The  network’s  training  phase  occurs  in  a  supervised  manner.  That  is,  a  training 
set  composed  of  input  vectors  v15 ...,  vN  and  their  associated  target  values  t  =  [tx  •  •  •  tN  ] 
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must  be  available.  To  begin,  the  center  vectors  are  placed  at  the  location  of  the  input 
vectors  themselves.  That  is,  n,  =  v;  for/  =  1,...,  N  .  The  NxN  matrix  of  hidden  layer 
neuron  activations  A  is  then  computed  in  which  the  elements  of  A  are 


(18) 


The  weight  vector w  =  [w,  •  •  •  wN~\  satisfies  the  equation 

Aw  =  t  (19) 

and  can  be  written  as 

w  =  ATt  (20) 

where  A'  represents  the  pseudo  inverse  of  A. 

Following  the  training  phase,  new  observations  can  be  presented  to  the  network  to 
generate  outputs.  Mathematically,  an  output  of  the  RBFNN  for  a  given  A-- feature  vector 
v  can  be  calculated  using 

y(v)  =  Yjwn  expf--^yXh  (21) 

n= 1  v  ^ n  k= 1  J 

where  wn  is  the  weight  of  the  nlh  center  vector,  vk  is  the  kth  feature  of  v,  and  //]"'  is  the  kth 
feature  of  the  nlb  center  vector. 

There  are  many  benefits  for  employing  ANNs  in  a  study.  In  fact,  the  power  of 
utilizing  multilayer  ANNs  comes  from  the  fact  that  any  continuous  function  can  be 
implemented  in  a  three-layer  neural  net  provided  that  there  is  a  sufficient  number  of 
hidden  layer  neurons  and  the  proper  nonlinear  activation  functions  are  chosen  [61]. 
Haykin  [63]  outlines  9  useful  properties  and  capabilities  for  ANNs: 
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1 .  The  network’s  nonlinear  nature  is  important  if  the  underlying  relationships 
between  the  components  of  the  input  signal  are  nonlinear. 

2.  The  ANN  utilizes  training  samples  to  leam  an  input-output  mapping  that  it 
applies  to  new  samples. 

3.  ANNs  are  flexible  to  changes  in  the  operating  environment. 

4.  As  it  relates  to  pattern  classification,  ANNs  offer  infonnation  about  which 
pattern  to  choose  as  well  as  a  level  of  confidence  in  making  that  choice. 

5.  The  network  manages  contextual  information  naturally  since  each  neuron  is 
potentially  affected  by  all  of  the  other  neurons  in  the  network. 

6.  ANNs  are  considered  fault  tolerant  in  that  when  a  network’s  performance 
degrades,  it  does  so  gracefully  instead  of  catastrophically. 

7.  A  neural  network  is  suitable  for  real-time  application  in  very -large-scale- 
integrated  (VLSI)  technology  that  necessitates  describing  complex  behavior 
in  a  hierarchical  manner. 

8.  Though  there  are  many  types  of  ANNs,  their  analysis  and  design  are  universal 
across  domains  and  their  commonalities  greatly  expand  the  ability  to  share 
theories. 

9.  ANN’s  correlation  to  the  brain  facilitates  an  expansion  in  the  areas  of  neural 
computing  as  well  as  neurobiology. 

Given  all  of  their  beneficial  properties  and  capabilities,  there  are  still  many 
concerns  regarding  the  use  of  ANNs.  First,  neural  networks  are  normally  used  to  make 
predictions  about  a  system  rather  than  to  build  models  or  develop  any  underlying 
knowledge  about  that  system  [9].  Second,  there  exists  a  risk  of  overfitting  the  data  when 
using  ANNs.  Neural  networks  can  provide  a  near-perfect  fit  to  historical  or  training  data, 
but  can  be  poor  at  predicting  new  data  [9,  35].  A  third  issue  when  working  with  ANNs 
lies  in  adjusting  the  network’s  complexity.  The  large  number  of  free  parameters,  or 
weights,  in  the  network  creates  difficulty  in  finding  a  balance  between  choosing  too  few 
and  too  many  neurons  to  achieve  the  best  generalization  of  the  phenomenon  of  interest 
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[61].  Finally,  although  ANNs  can  be  a  powerful  and  fast  tool,  they  should  be  used  as  a 
supplement,  not  a  substitute,  to  standard  regression  and  designed  experiments  statistical 
tools  since  they  do  not  allow  fundamental  insights  into  the  underlying  system  mechanism 
that  produced  the  data.  Neural  networks  cannot  provide  the  solution  on  their  own  and 
should  be  integrated  into  a  consistent  system  engineering  approach  [9,  35,  61,  63]. 

2.4  Optimization  Approaches  Using  the  Mean  and  Variance  Models 

Given  that  there  are  available  mean  and  variance  models,  ju  and  cr2  respectively,  a 
number  of  dual  response  optimization  approaches  can  be  taken  to  identify  a  system’s 
robust  control  factor  setting.  Vining  and  Myers’  [6]  determined  robust  operating 
conditions  by  optimizing  the  primary  response  model  subject  to  a  constraint  on  the 
secondary  response  model.  In  a  smaller-the-better  case,  cr  is  restricted  at  some  specified 
value  a\  while  ju  is  minimized.  A  larger-the-better  situation  maximizes  ju  while 

controlling  a2  at  some  specified  value  aj  .  Finally,  in  a  target-is-best  case,  the  concern  is 
maintaining  ju  at  some  specified  value  /uT  while  a2  is  minimized.  These  three  optimization 
problems  are  shown  in  Table  2. 


Table  2.  Vining  and  Myers'  Dual  Response  Optimization  Scenarios 


Smaller-the-Better 

Larger-the-Better 

Target-is-Best 

Minimize  ju 

Maximize  ju 

Minimize  a2 

Subject  to  a 2  =  a 2 

Subject  to  cf  =  cr2 

Subject  to  ju  =  jUT 
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Lin  and  Tu  [46]  contended  that  Vining  and  Myers’  use  of  equality  constraints 
most  likely  eliminates  finding  better  global  solutions.  Lin  and  Tu  make  the  case  that  a 
better  solution  can  be  found  by  allowing  some  deviation,  or  bias,  of  the  mean  around  the 
target  value  jUT  while  keeping  the  variance  small.  Their  method  solves  for  optimal  control 

factor  settings  by  estimating  the  MSE  using  a  function  of  the  process  mean  and  variance. 
Their  MSE  criteria  for  the  three  experimental  scenarios  are  formulated  in  Table  3. 
Optimal  control  factor  settings  can  be  found  by  solving  the  optimization  problem 

Minimize  MSE 

(22) 

Subject  to  x  e  D 

where  D  is  the  experimental  design  space. 


Table  3.  Lin  and  Tu's  Dual  Response  Optimization  Scenarios 


Smaller-the-Better 

Larger-the-Better 

Target-is-Best 

MSE  =  <j2+ju2 

MSE  =  &2  -ft2 

MSE  =  a2  +  ( jU  -  jury 

One  criticism  of  Lin  and  Tu’s  target-is-best  method  is  that  there  is  no  restriction 
on  how  far  the  mean  process  response  may  deviate  from  the  target  value  and,  as  a  result, 
may  be  deficient  if  it  is  critical  to  maintain  the  mean  close  to  the  target  [65].  In  situations 
such  as  these,  Copeland  and  Nelson  [66]  recommended  obtaining  a  solution  in  which  //  is 

within  a  specified  distance  (A  )  of  /uT  .  They  endorsed  minimizing  a  subject  to 

[ju-  jUT)2  <  A2 .  The  additional  smaller-the-better  and  larger-the-better  instances  are 

shown  in  Table  4  where  A  ,  is  a  maximum  allowable  value  for  d2 . 

(7 
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Table  4.  Copeland  and  Nelson's  Dual  Response  Optimization  Scenarios 


Smaller-the-Better 

Larger-the-Better 

Target-is-Best 

Minimize  ju 

Minimize  -  ju 

Minimize  d2 

Subject  to  d2  <  A  2 

Subject  to  d2  <  A^2 

Subject  to  (ju- jUt)~  <  A 2 

Ding  et  al.  [67]  suggested  a  weighted  MSE  ( WMSE )  approach  by  utilizing  the 
convex  combination  of  the  mean  and  variance  functions.  Their  proposal  minimizes  the 
WMSE s  in  Table  5  where  A  e  [0,1] . 


Table  5.  Ding  et  al.’s  Dual  Response  Optimization  Scenarios 


Smaller-the-Better 

Larger-the-Better 

Target-is-Best 

WMSE  =  Aju2  +  (\-X)d 2 

WMSE  =  -Ajj2  +(\-A)d2 

WMSE  =  2  (ju  -  jur  f  +  (1  -  A)  d2 

The  methodologies  of  Vining  and  Myers,  Lin  and  Tu,  Copeland  and  Nelson,  and 
Ding  et  al.  are  useful  when  a  single  optimal  solution  is  necessary.  Koksoy  and 
Doganaksoy  [68]  present  a  flexible  nonlinear  multi-objective  approach  by  considering  the 
secondary  response  as  another  primary  response.  They  claim  that  the  restriction  placed 
upon  the  secondary  response  may  exclude  better  conditions.  Their  method,  which  only 
focuses  on  the  smaller-the-better  and  larger-the-better  problem  structures,  allows  further 
insight  into  the  RPD  problem  by  exploring  trade-offs  between  the  mean  and  variance 
responses.  These  formulations  are  shown  in  Table  6.  Solving  these  problems  results  in 
finding  a  string  of  Pareto  alternative  solutions  in  some  region  of  interest  R  that  jointly 
optimize  ju  and  cr 2 .  That  is,  it  is  impossible  to  improve  ju  without  making  cr2  worse  and 
vice  versa. 


24 


Table  6.  Koksoy  and  Doganaksoy’s  Dual  Response  Optimization  Scenarios 


Smaller-the-Better 

Larger-the-Better 

|  Minimize  ju,  Minimize  <j2j 

jMaximize  ju,  Minimize  cr  j 

Subject  to  x  €  R 

Subject  to  x  €  R 

2.5  Multi-Response  Robust  Parameter  Design  Approaches 

In  the  multi-response  RPD  problem,  the  objective  is  to  find  the  optimal  control 
parameter  levels  that  return  average  responses  close  to  their  target  values  while 
minimizing  the  variance  of  each  response.  The  literature  offers  several  methods  for 
optimizing  multi-response  problems.  These  methods  involve  the  use  of  desirability 
functions  [14-18],  loss  functions  [19-23],  principal  component  analysis  (PCA)  [24-27], 
distance  metrics  [28,  29],  and  MSE  criterion  [30-32].  A  majority  of  these  techniques 
transform  the  quality  characteristics  into  new  response  variables  in  order  to  reduce  the 
dimension  of  the  optimization  problem.  Typically,  the  transformations  convert  the  output 
responses  to  a  single  response. 


2.5.1  Desirability  Functions 

Derringer  and  Suich  [14]  adopted  Harrington’s  [69]  use  of  desirability  functions 
which  map  each  of  the  estimated  responses  yt  into  a  desirability  value  cL  where  0  <  cL  <  1 . 

An  overall  system  desirability  D  is  then  generated  by  combining  the  J  individual 
desirabilities  via  the  geometric  mean  as  in  Equation  (20). 

(  J  ^ 


D  = 


(23) 


V  7=i  J 
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An  optimal  operating  point  for  the  set  of  responses  is  then  found  by  maximizing  D. 
Given  the  actual  desire  for  the  individual  output  responses  in  the  experiment,  such  as  a 
minimum  or  maximum  value,  the  desirability  values  can  be  defined  as  in  Table  7. 


Table  7.  Derringer  and  Suich's  Desirability  Functions 


Smaller-the-Better 


Larger-the-Better 


Target-is-Best 


d,  = 


1  yt  <  L 

/  n  r 

U,-y, 


L  <  y.  <  U  d, 
S.V-L,)  '  '  ' 

0  y,>U, 


y,  <  A 


A -A 


A  d  = 

^  i  J 

1  yi>Ui 


fy-  A^ 

vA-Ay 

ru-y' 

kU-t.j 


A  ^  yt  ^  A 

A  <y,*Ui 

y.  <  Lt 

or 

9i>Ui 


The  minimum  and  maximum  allowable  values  for  y(  are  denoted  A  and  A, , 
respectively.  These  points  can  also  represent  levels  at  which  permitting  either  j).  <  L,  or 
yt  >  Uj  adds  very  little  value  to  the  overall  process.  Also,  A  is  the  desired  target  value  of 
j>,.  between  A  and  A; .  The  exponents  r  and  s  operate  as  shape  parameters  for  the 
desirability  function.  A  large  value  for  r  or  s  puts  greater  importance  on  the  response 
values  being  closer  to  the  respective  target.  Smaller  values  imply  that  the  desirability 
value  is  large  even  if  the  response  is  far  from  its  target  value.  The  desirability  functions 
for  each  of  Taguchi’s  three  experimental  cases  are  illustrated  in  Figure  3.  Finally,  to  find 
the  optimal  process  settings,  D  is  maximized  with  respect  to  the  controllable  factors. 
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Smaller-the-Better  Larger-the-Better  Target-1  s-Best 


Figure  3.  Desirability  Functions  for  Taguchi’s  Three  Experimental  Cases 


There  are  several  variations  to  Derringer  and  Suich’s  methodology.  Instead  of  the 
geometric  mean,  Park  [15]  advocated  the  use  of  the  harmonic  mean  of  the  /  desirabilities. 
Del  Castillo  et  al.  [16]  presented  modified  desirability  functions  that  are  everywhere 
differentiable.  Park  and  Park  [17]  introduced  a  weighted  desirability  function  approach 
that  allows  for  varying  degrees  of  importance  to  be  applied  to  the  different  responses. 

The  weighted  desirability  formulation  is 

(  j  \Yj 

dw=  IK'  (24) 

V  J=l 


where  w .  >  0  is  the  weight  of  the /h  response  and  ^  w .  =  J  .  Wang  et  al.  [  1 8]  proposed  a 

j= i 


robust  desirability  function. 


27 


2.5.2  Loss  Functions 


Pignatiello  [19]  focused  on  Taguchi’s  target-is-best  experimental  scenario  and 
based  his  method  on  the  multivariate  quadratic  loss  function 

^[y(x)]  =  (y(x)-T)'c(y(x)-T)  (25) 

where  y(x)  is  a  vector  of  responses  for  parameter  setting  x,  t  is  the  target  response  vector, 
and  C  is  a  cost  matrix.  The  optimal  parameter  setting  is  then  be  found  by  minimizing  the 
expected  loss  defined  by 

£[My(x)]]  =  (£[y(x)]-T)  C (F [y (x)]  - t)  +  trace (C  •  Cov [y(x)])  (26) 

where  E [y(x)]  and  Cov[y(x)]  are  the  respective  mean  vector  and  covariance  matrix  of 

the  responses  at  parameter  setting  x.  Many  authors  have  modified  Pignatiello’s  approach. 
Ames  et  al.  [20]  developed  loss  functions  focusing  on  the  individual  responses  being  on 
target,  but  having  no  consideration  for  the  correlation  structure  between  the  responses. 
Vining  [21]  modified  Equation  (26)  to  use  the  mean  and  variance-covariance  structure  of 
the  predicted  responses  y(x)  instead  of  the  mean  and  variance-covariance  structure  of  the 
actual  responses  y(x) .  This  approach  considered  the  prediction  quality  as  well  as  the 
correlation  structure  of  the  responses.  Romano  et  al.  [22]  adopted  Vining’s  multi¬ 
response  quality  loss  function  and  minimized  the  expected  total  loss  subject  to  lower  and 
upper  bound  constraints  placed  on  the  means  and  variances  of  the  individual  responses. 
Ko  et  al.  [23]  integrated  the  strengths  of  Pignatiello  and  Vining’s  approaches  to  minimize 
the  expected  future  loss.  One  drawback  to  the  loss  function  approach  is  that  it  may  be 
difficult  to  define  the  cost  matrix  C  [70]. 
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2.5.3  Principal  Component  Analysis 

Su  and  Tong  [24]  grounded  their  strategy  on  transfonning  the  nonnalized  quality 
losses  into  a  set  of  uncorrelated  components  via  PCA.  Salmasnia  et  al.  [25]  first 
transformed  principal  component  models  into  a  desirability  function.  They  then  found  an 
optimal  solution  by  maximizing  the  overall  desirability  of  the  selected  principal 
components  within  the  desired  region  of  the  normalized  means  and  standard  deviations  of 
the  original  responses.  Paiva  et  al.  [26]  combined  the  PCA  and  MSE  approaches  into  a 
Multivariate  Mean  Square  Error  (MMSE)  measure.  Gomes  et  al.  [27]  expanded  on  Paiva 
et  al.’s  approach  and  presented  the  Weighted  Multivariate  Mean  Square  Error  (WMMSE) 
to  appropriately  weight  the  individual  responses  in  the  MMSE  approach. 

2.5.4  Distance  Metrics 

Khuri  and  Conlon  [28]  considered  the  Mahalanobis  distance  between  a  vector  of 
each  response  function  and  a  corresponding  vector  of  their  optimum  function  values.  The 
robust  solution  across  the  set  of  responses  is  the  x*  in  the  experimental  region  that 
minimizes  this  distance.  Govindaluri  and  Cho  [29]  decoupled  the  J  individual  response 
MSE  functions  from  the  expected  total  loss.  They  then  found  the  setting  that  minimized 
the  distance  between  the  vector  of  individual  response  MSE  functions  and  the  vector  of 
ideal  MSE  values.  Chiao  and  Hamada  [71]  modeled  the  mean  and  covariance  structure 
of  the  assumed  multivariate  normal  responses  and  then  maximized  a  “proportion  of 
conformance”  measure  defined  as  the  probability  that  the  /  responses  jointly  meet  their 
respective  specification  limits. 
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2.5.5  Mean  Square  Error  Criterion 

Koksoy  and  Yalcinoz  [30],  Koksoy  [31],  and  Koksoy  [32]  extended  Lin  and  Tu’s 
[46]  MSE  criterion  to  the  multi-response  robust  design  case.  Koksoy  and  Yalcinoz  [30] 
promoted  the  minimization  of  the  weighted  summation  of  the  individual  MSE  functions 


Minimize  y^  WMSE i 

i= i 

Subject  to  xeR 


(27) 


where  Wj  is  the  weight  of  the  /h  MSE  function,  ^  If 7  =  1 ,  and  R  is  the  region  of  interest. 


7=1 


Koksoy  [31]  recommended  solving  the  following  multi-objective  optimization  problem: 

(28) 


Minimize  { MSE] , MSE2 MSEj } 
Subject  to  xeR 


In  the  non-trivial  multi-objective  optimization  problem,  a  single  solution  does  not  exist 
that  simultaneously  optimizes  each  objective.  Therefore,  a  list  of  Pareto  optimal 
solutions  in  which  an  improvement  to  one  objective  causes  degradation  to  at  least  one 
other  objective  is  generated  for  the  decision  maker.  All  Pareto  solutions  are  considered 
equally  good.  Koksoy  [32]  further  proposed  solving  the  optimization  problem 
Minimize  MSE .  for  1  <  /'  <  J 


Subject  to  MSEj  =  MSEj0  for  i  =  1,2,...,  J;  i  ^  j 
X  <E  R 


(29) 


where  MSE,,,  are  specified  values  for  each  MSE  function.  By  successively  changing  the 
specified  constraint  values,  a  string  of  solutions  is  generated  rather  than  a  single  solution. 
This  allows  for  an  improved  understanding  of  the  problem  by  examining  the  trade-offs 
that  must  be  considered  to  obtain  a  compromised  solution. 


30 


2.6  The  Delta  Method 


Situations  surface  where  interest  lies  in  the  distribution  of  some  nonlinear 
function  of  a  random  variable  and  not  necessarily  the  distribution  of  the  random  variable 
itself.  With  that,  the  concern  then  turns  to  the  properties  of  the  function  of  the  random 
variable.  In  particular,  how  can  the  variance  of  the  function  of  the  random  variable  be 
estimated?  One  such  technique,  known  as  the  Delta  Method,  utilizes  a  Taylor  series 
approximation  to  obtain  reasonable  estimates  for  the  mean  and  variance  of  the  function  of 
a  random  variable. 

Though  the  original  author  of  the  Delta  Method  is  unknown,  an  article  by  Robert 
Dorfman  in  the  1938  journal  Biometric  Bulletin  is  credited  as  the  earliest  use  of  the  “S- 
method.”  According  to  Ver  Hoef  [72],  Dorfman  proposed  the  technique  to  approximate 
the  variance  of  a  nonlinear  function  of  multiple  random  variables.  Ver  Hoef  also 
reproduced  Dorfman’s  contribution  in  which  he  comments  that,  if/is  a  linear  function, 
then  the  (5-method  is  exact.  He  also  states  that,  if/“does  not  deviate  sharply  from  linear,” 
the  ^-method  gives  a  good  approximation. 

Since  its  origin,  scientists  and  statisticians  across  a  variety  of  fields  have  utilized  a 
version  of  Dorfman’s  original  (5-method.  Chapra  and  Di  Toro  [73]  extended  the  Delta 
Method  to  modeling  water  quality  and  estimating  stream  reaeration,  production,  and 
respiration  rates.  Durbin  et  al.  [74]  utilized  the  Delta  Method  to  derive  a  transformation 
of  non-nonnally  distributed  DNA  microarray  data  to  stabilize  the  asymptotic  variance 
over  the  full  range  of  data.  Powell  [75]  focused  on  providing  variance  approximations 
for  common  parameters  used  by  avian  ecologists  such  as  annual  bird  population  growth 
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and  mean  annual  density  of  bird  species.  White  [76]  has  also  developed  a  Windows- 
based  software  called  MARK  which  utilizes  the  Delta  Method  to  aid  in  the  parameter 
estimation  of  marked  animals  when  they  are  re-encountered  via  dead  recoveries,  live 
recaptures,  or  radio  tracking. 


2.6.1  Univariate  Case 

Recall  from  calculus  [77]  that  if  a  function  /:  R  -»  R  has  derivatives  of  order  n, 


that  is  f°'\x)  = 


dnf{x) 

dxn 


exists,  then  the  Taylor  series  expansion  of f  centered  at  some 


constant  a  is  defined  as 


1,-0  n\ 


Consequently,  the  second-order  Taylor  series  approximation  of /  centered  at  a  is 


(30) 


fix)  ~  f{a)  +  f'(a)(x  -  a)  +  \f\a){x  -a)~.  (31) 

Now  consider  Y  =  f  ( X )  as  a  function  of  the  normally  distributed  random  variable 
X.  Casella  and  Berger  [78]  consider  estimating  the  mean  and  variance  of  Y  when  the 
mean  and  variance  of  X are  known  parameters.  That  is,  the  interest  is  in  finding  E  [T]  and 

Var[Y ]  given  that E[X]  -  fix  and  Var(X)  =  a\.  Following  from  Equation  (31),  the 


second-order  Taylor  series  approximation  of  Y  centered  at  the  point  a  =  jux  is  defined  as 


7  =  f(X)  *  f(<ux)  +  n<ux)(X- Mx)  +  ^f"(Mx)(X-Mx)2 . 


(32) 


Therefore,  applying  the  expectation  operator  to  Equation  (32)  generates  an  estimate  for 
the  mean  of  Y: 
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(33) 


Similarly,  applying  the  variance  operator  to  Equation  (32)  generates  an  estimate  for  the 
variance  of  Y : 

Var[Y]«f^xfa2x  t}  (34) 

Derivations  of  Equations  (33)  and  (34)  are  shown  in  Appendix  A. 

2.6.2  Multivariate  Case 

In  the  multivariate  case  [79],  let  S'  be  a  nonempty  set  in  Mn  and  let  f:  S  ->  E.  For 
x  e  S  ,  if  the  gradient  vector  V/(x)  and  hessian  matrix  H(x)  exist,  then  the  multivariate 
second-order  Taylor  series  approximation  of/ centered  at  some  constant  vector  a  is 

/ (x)  *  /(a)  +  V/  (a)'(x  -  a)  +  \  (x  -  a)'H(a)(x  -  a) .  (35) 

Now  consider  a  function  /:  Rn  E  and  a  n-dimensional  normally  distributed 
random  vector  X  with  mean  vector  px  and  covariance  matrix  Zx  .  The  second-order 

Taylor  series  approximation  of  the  function  Y  =  f  (X)  centered  at  the  vector  a  =  px  is 

/(x)«/o»x)+v/oixy(x-|ix)+i(x-|ixyH(|ixxx-|ix).  (36) 

The  estimated  mean  and  variance  of  Y  are  then  defined  respectively  as 

£[Jl*/(M  +  MH^x)2x)  (37) 

and 

Far[7]  =  V/(^)'i;xV/(fc)+}fr(H(fc)ExHK)Ix).  (38) 

Derivations  of  Equations  (37)  and  (38)  are  shown  in  the  Appendix  A. 
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III.  Extending  the  Combined- Array  Response  Surface  Model  Approach 


3.1  Introduction 

The  first  objective  of  this  dissertation  is  to  extend  the  combined-array  RSM 
approach  that  relies  exclusively  on  the  low-order  polynomial  models  discussed  in  Section 
2.3.1.  Since  more  accurate  predictive  response  surface  models  result  in  better  RPD 
solutions  [33],  a  methodology  will  be  developed  that  utilizes  the  non-linear  Kriging  and 
RBFNN  models  in  place  of  the  polynomial  models.  From  there,  the  mean  and  variance 
of  a  second-order  Taylor  series  approximation  of  the  Kriging  and  RBFNN  models  will  be 
calculated  via  the  Multivariate  Delta  Method.  Finally,  an  existing  optimization  problem 
that  employs  these  approximations  will  be  solved  to  identify  the  robust  control  parameter 
setting.  Henceforth,  this  procedure  is  referred  to  as  the  combined-array  Multivariate 
Delta  Method  approach,  or  simply  MDM. 

The  rest  of  Chapter  III  is  organized  as  follows.  Section  3.2  outlines  the  proposed 
MDM  methodology.  Section  3.3  uses  two  case  studies  to  compare  the  combined-array 
MDM  approach  to  the  combined-array  RSM  approach  and  the  combined-array  stochastic 
emulator  approach.  Finally,  Section  3.4  summarizes  this  chapter. 

3.2  The  Combined-Array  Multivariate  Delta  Method  (MDM)  Approach 

This  section  outlines  the  methodology  behind  the  proposed  MDM  approach.  In 
Section  3.2.1,  the  mean  and  variance  models  for  the  second-order  Taylor  series 
approximations  to  the  Kriging  and  RBFNN  models  are  developed.  In  Section  3.2.2,  the 
combined-array  MDM  approach  is  outlined  against  the  combined-array  RSM  approach 
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and  the  combined-array  emulator  approach.  A  cross-validation  procedure  for 
determining  an  appropriate  network  structure  for  the  RBFNN  is  highlighted  in  Section 
3.2.3. 


3.2.1  Mean  and  Variance  Models  via  Taylor  Series  Approximation 

The  MDM  approach  uses  the  same  distributional  assumptions  regarding  the  noise 
variables  and  the  model’s  error  as  the  RSM  approach  does.  To  recap,  the  model’s  error  e 
is  normally  distributed  with  a  zero  mean  and  a  constant  variance  cr?  .  Also,  it  is  assumed 

that  the  noise  factors  are  uncorrelated  and  that  z  ~  A(pz ,  Zz ) . 


First,  let  j)(v) ,  where  v  = 


x 


represent  the  Kriging  or  RBFNN  model  of  the 


LZJ 

system.  The  second-order  Taylor  series  approximation  o f y ( v )  about  the  vector 
fxl 


a  =  Hv 


is  defined  as 


r2(i)(V))  =  i)(Hv)  +  Vi)(llvy(V-l1v)  +  KV_llvyH(Pv)(v-llv)  +  'e  (39) 
where  Vj)(pv)  and  H(pv )  are  the  gradient  vector  and  Hessian  matrix  of  j)(v)  evaluated  at 
pv .  By  way  of  the  Multivariate  Delta  Method  [78],  the  estimated  mean  and  variance  of 
T2  (y(v))  are  then  calculated  as 

E[T2  (  v(v))]  =  Kll,)  +  l‘r(  H(|lv)S,)  (40) 

and 

Var[_T2  (f(v))]  =  Vy(pv)'Zv  Vv(pv)  +  yP-(H(pv)EvH(pv)Zv)  +  cr^  ,  (41) 
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Since  the  Kriging  model  predicts  the  exact  response  for  observed  training  vectors,  its 
MSE  is  equal  to  zero.  The  MSE  of  the  RBFNN  model  is  discussed  in  Section  3.2.3. 

Finally,  in  a  manner  similar  to  the  expressions  for  the  mean  and  variance  of  the 
standard  and  extended  quadratic  models,  Equations  (42)  and  (43)  are  in  terms  of  only  the 
control  factors  x  and,  as  such,  they  can  be  used  to  approximate  the  mean  and  variance  of 
the  system  anywhere  in  the  control  design  space.  Given  these  mean  and  variance 
estimates,  a  dual  response  optimization  approach  is  then  solved  to  locate  a  robust  control 
parameter  setting.  A  subset  of  the  optimization  formulations  is  discussed  in  Section  2.4. 

Below,  Sections  3. 2. 1.1  and  3.2. 1.2  derive  the  gradient  vector  and  Hessian  matrix 
of  the  Kriging  and  RBFNN  metamodels. 
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3.2. 1.1  Gradient  Vector  and  Hessian  Matrix  of  the  Kriging  Model 

Let  K  =  rx  +  rz .  The  kth  element  of  the  gradient  vector  for  the  Kriging  model 


output  in  Equation  (15)  for  p  =  2  is  defined  as 


^  =  ^>R-‘(y-l/?)  for k  = 


dv,  dv 


(44) 


where 


<9r(v) 

dv. 


=  -20,. 


(v,  -v!'>)exp 


-Z«.( 


v  —  v(1)) 

m  m  / 


{yk  -  V\N)  )  eXP  “Z  °m  (■ Vn,  -  vr  )2 


m= 1 


is  the  first-order  partial  derivative  ofr(v)  with  respect  to  input  feature  k. 

The  ijth  element  of  the  Hessian  matrix  for  the  Kriging  output  is  defined  as 
d2y  32r'(v) . 


dvfvj  dvfivj 


-R_1(y-1/?)  for  i,j  = 


where 


(45) 


(46) 


S2r(v) 

dvfivj 


20 


(2 0t  (v;.  -  v®  )2  - 1)  exp  ~Z  0m  (vm  -  v™ ) 


m= 1 


(2  0,  (v,.  -v,(A,))2  -l)exp  ~Yd0m  (vm -v\r) 


(v,.  -  v'1'  )(v j  -  vj1’ )  exp  -Z  0m  (vm  -  v® ) 

m= 1 


4^, 


for  i  =  j 


(y.-v^Jfv.-v^Jexp 


-TO  (v  -vw)2 

m  \  m  m  ) 


n= 1 


for  i  j 


(47) 


is  the  second-order  partial  derivative  ofr(v)  with  respect  to  input  features  i  and  j. 
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3.2. 1.2  Gradient  Vector  and  Hessian  Matrix  of  the  RBFNN  Model 
The  kth  element  of  the  gradient  vector  for  the  RBFNN  output  in  Equation  (21)  is 
defined  as 


dy  _ 

dvk~  h  o>H 

The  ifh  element  (for/,/  =  1,. 
defined  as 


exp 


Z<Jn  '"=1 


fork  =  1  , 


(48) 


,.,K  )  of  the  Hessian  matrix  for  the  RBFNN  output  is 


f2v 

dv:dv,. 


n= 1 

N  W, 


<J„ 


exp 


1 

m= 1 


2a; 


n= i 


cr„ 


exp 


~Tl(v.  -M™)2 

Z<-^n  m=\ 


for  /  =  j 
for  /  ±  j 


(49) 


3.2.2  RSM  Approach  vs.  Emulator  Approach  vs.  MDM  Approach 

The  main  steps  for  the  RSM  approach,  the  emulator  approach,  and  the  MDM 
approach  are  outlined  in  Figure  4.  This  figure  also  highlights  the  experimental  designs 
(DOE)  and  the  modeling  efforts  associated  with  each  approach.  The  RSM  and  MDM 
approaches  each  require  the  development  of  one  DOE  and  one  modeling  effort.  The 
emulator  approach,  on  the  other  hand,  necessitates  the  generation  of  two  DOEs  and  three 
individual  modeling  efforts. 

In  order  to  compare  the  results  across  the  three  different  approaches,  each  j)(x,  z) 
was  built  from  a  combined-array  DOE  in  which  the  coded  ±1  levels  of  each  noise  factor 
z;  were  set  to  //,  ±  3cr  .  This  ensures  that  the  emulator  remains  valid  during  the  sampling 
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phase  [13].  Consequently,  setting  the  levels  of  the  noise  factors  in  this  manner  implies 
thatZz  =  (y)2 1, 

The  last  step  of  each  approach  is  to  identify  the  robust  control  setting  by  utilizing 
the  appropriate  mean  and  variance  models  within  an  optimization  problem.  Due  to  its 
simplicity  and  the  manner  in  which  it  balances  being  on  target  with  minimal  variance, 

Lin  and  Tu’s  MSE  approach  in  Table  3  was  used  in  this  research.  The  optimal  settings  x* 
were  found  using  MATLAB’s  fmincon  function  was  used  within  a  greedy  randomized 
adaptive  search  procedure  (GRASP)  heuristic  that  solves  the  problem  from  a  number  of 
different  starting  points  and  chooses  the  best  overall  solution  [80].  This  procedure  can  be 
slow  because  it  seeks  to  avoid  getting  trapped  at  a  local  minimum.  This  is  a  concern 
since  fmincon  does  not  guarantee  convergence  to  a  global  minimum  for  a  potentially 
nonconvex  problem  [81]. 
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Figure  4.  Combined-Array  RPD  Approaches:  RSM  vs.  Emulator  vs.  MDM 


3.2.3  Choosing  the  RBFNN  Structure 

Although  RBFNNs  have  many  advantages,  one  disadvantage  they  have  is  that,  as 
the  size  of  the  training  set  increases,  the  number  of  hidden  neurons  that  are  needed  also 
increases.  This  subsequently  increases  the  time  necessary  to  train  the  network.  Another 
drawback  to  their  use  involves  the  selection  of  the  network’s  architecture  [61].  If  too 
many  neurons  are  used,  then  the  overall  generalization  of  the  network  will  be  deficient. 
On  the  other  hand,  the  network  will  not  be  able  to  sufficiently  learn  the  training  data  if 
too  few  neurons  are  selected.  Therefore,  a  cross-validation  (CV)  procedure  was  used  to 
determine  the  structure  of  the  RBFNN  so  that  the  resulting  function  was  well-generalized 
with  minimal  risk  of  over-fitting  the  experimental  data.  The  structural  parameters  of 
concern  here  were  the  number  of  hidden  layer  neurons  used  in  the  network  and  the  spread 
parameter  a  of  the  hidden  layer  neurons. 

Given  an  experimental  design  with  N  runs,  a  set  of  mn  hidden  layer  sizes  and  a  set 
ofmCT  spread  values  were  first  defined.  Then,  M  =  mn  x  ma  combinations  of  RBFNN 
structure  pairs  were  generated  such  that  each  pairing  was  composed  of  a  hidden  layer  size 
(n)  and  a  spread  value  (cr).  Next,  a  CV  procedure  was  performed  across  the  set  of 
structures  and  the  MSE  was  recorded  for  each  (n,  a)  pair.  The  size  of  the  experimental 
design — this  research  used  designs  of  25,  81,  and  256  runs — detennined  which  CV 
method  was  used.  For  N <  25,  leave-one-out  CV  was  used.  For  computational 
efficiency,  ten  rounds  of  10-fold  CV  were  used  for  N  >  25.  The  final  neural  network  was 
trained  on  the  complete  design  and  structured  via  the  (n,  a)  pairing  that  produced  the 
minimum  MSE  (or  minimum  average  MSE  for  larger  datasets).  For  the  examples  used  in 
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this  research,  this  procedure  led  to  well-generalized  networks  with  a  reduced  risk  of  over¬ 
fitting  the  data  since  only  20-50%  of  the  total  number  of  training  vectors  were  routinely 
chosen  as  neuron  centers.  The  chosen  spread  parameter  ranged  from  2  to  10.  Since  the 
design  space  was  limited  between  -1  and  1,  these  values  for  the  spread  parameter  led  to 
smooth  functions.  The  RBFNNs  were  trained  using  MATLAB’s  newrb  function. 

The  MSE  of  the  RBFNN  that  is  used  to  estimate  cr]  in  Equation  (43)  is 

/(H-")  (50) 

;=1  / 

where  e,  is  the  prediction  error  of  the  /lh  design  point.  The  value  n,  which  is  the  number  of 

neurons  in  the  trained  network’s  hidden  layer,  corresponds  to  the  number  of  estimated 
weight  parameters  in  the  RBFNN. 

3.3  Application  and  Results 

In  this  section,  the  MDM  approach  is  applied  to  two  case  studies.  Section  3.3.1 
provides  a  proof-of-concept  demonstration  of  the  MDM  approach  using  a  synthetic  case 
study.  Section  3.3.2  applies  the  MDM  approach  to  a  computer  simulation. 

A  popular  method  for  comparing  different  analysis  techniques  for  simulation 
studies  is  through  the  use  of  the  Root  Mean  Squared  Error  ( RMSE )  [43,  5 1,  53,  54,  60]. 
Since  the  analysis  is  of  known  models,  the  model  predictions  ym  can  be  compared  to  the 

known  true  values  ym  for  a  test  set  ofM  =  200  random  validation  points  via 

i  M 

MSE-^y.-yJ  .  (51) 
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A  model’s  RMSE  can  then  be  compared  to  the  range  of  the  true  responses  to  gain  insight 
into  its  accuracy  [11]. 

3.3.1  A  Synthetic  Case  Study 

In  this  case  study,  data  was  generated  from  the  truth  model  in  Equation  (5 1).  By 
knowing  the  truth  model,  the  system’s  true  mean  and  variance — shown  in  Equations  (52) 
and  (53) — can  be  used  as  reference  points  for  our  modeling  efforts.  The  truth  model  is 
employed  with  the  control  factory  e  [-1,1]  and  the  noise  factor  z,  ~A(/r  =  0,cr  =  l).  The 

goal  of  this  RPD  study  was  to  locate  the  appropriate  operating  point  that  resulted  in  a 
mean  response  of  8  with  a  minimal  variance.  A  5  full  factorial  design  was  used  to 
generate  experimental  data. 


=  9//(l  +  exp(l-5x1))  +  3z12  - 4x,2z, 

(52) 

£  [v]  =  9/(  1  +  exp  (1  -  5xj ))  +  3 

(53) 

Var  [  v]  =  16xj4  +  18 

(54) 

The  MDM  approach  using  Kriging  (KR)  and  RBFNN  (RBF)  models  was 
compared  to  the  RSM  approach  that  uses  the  four  polynomial  models  in  Section  2.3.1: 
y std  (x’z)’  y nn  (x>  z)  ,  yCNN  (x> z)  >  and  yCcN  (x> z)  •  Henceforth,  these  models  will  simply  be 
referred  to  as  Std,  NN,  CNN,  and  CCN. 

Response  surface  graphs  for  the  six  modeling  efforts  are  shown  alongside  the 
truth  model  in  Figure  5.  Visual  inspection  shows  that  the  Std,  NN,  and  CNN  models  only 
provide  the  general  trend  of  the  data.  This  was  expected  due  to  the  highly  non-linear 
nature  of  Equation  (52).  The  CCN  model  offers  a  markedly  improved  representation  of 
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the  non-linearity  of  the  system  over  its  predecessors.  Finally,  the  KR  and  RBF  models 
represent  the  true  input-output  relationship  of  the  system  very  well.  This  information  is 
summarized  in  Table  8. 


Figure  5.  Response  Surface  Comparison  for  the  Synthetic  Case  Study 


Table  8.  Quality  of  Models  for  the  Synthetic  Case  Study 


Model 

RMSE 

%  of  Range 

Std 

1.62 

11.39% 

NN 

1.13 

7.97% 

CNN 

1.13 

7.97% 

CCN 

0.79 

5.58% 

RBF 

0.31 

2.18% 

KR 

0.19 

1.37% 

The  mean,  variance,  and  MSE  models  that  correspond  to  the  Truth,  Std,  NN, 


CNN,  CCN,  KR,  and  RBF  response  models  are  depicted  in  Figure  6.  The  locations  of  the 


estimated  robust  points  are  also  plotted.  In  this  case,  the  NN,  CNN,  and  CCN  mean 
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models  are  nearly  identical  and  are  plotted  on  top  of  one  another.  Similarly,  the  NN  and 
CNN  variance  models  are  graphically  the  same.  Of  the  six  modeling  efforts,  the  KR  and 
RBF  mean  models  are  very  good  representations  of  the  true  mean  response.  Also,  the 
CCN,  KR,  and  RBF  variance  models  provide  the  best  depictions  of  the  system’s  true 
variance.  In  fact,  the  RBF  variance  model  is  nearly  identical  to  the  system’s  true  variance 
model.  The  Std,  NN,  and  CNN  mean  and  variance  models  are  poor  relative  to  their 
competitors. 


Mean  Model 


Variance  Model 


MSE  Model 


- True  Model  - NN  Model 

True  Robust  Point  -+  NN  Robust  Point 

-  Std  Model  - CNN  Model 

Std  Robust  Point  O  CNN  Robust  Point 


CCN  Model  - RBF  Model 

□  CCN  Robust  Point  *  RBF  Robust  Point 

- KR  Model 

0  KR  Robust  Point 


Figure  6.  Mean,  Variance,  and  MSE  Models  for  the  Synthetic  Case  Study 


Table  9  summarizes  the  predicted  and  realized  means,  variances,  and  MSEs  at 
each  robust  point.  Confirmation  experiments  were  performed  in  which  2,000  Monte 
Carlo  simulations  of  the  truth  model  were  run  at  the  estimated  robust  points.  Common 
random  numbers  were  used  to  allow  an  apples-to-apples  comparison.  It  can  be  observed 
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that,  not  only  is  the  RBF  robust  point  the  closest  to  the  true  robust  point,  but  its 
predictions  are  indicative  of  the  system’s  actual  performance.  An  interesting  note  is  that, 
without  knowledge  of  the  system’s  truth  model — which  is  typically  never  known  for  a 
simulation — and  based  solely  on  the  predictions  provided  in  Table  9,  initial  conclusions 
would  have  chosen  the  Std  robust  point  as  the  best  operating  point  since  the  prediction  of 
the  system’s  mean  response  is  exactly  on  target  and  the  prediction  of  the  system’s 
variance  is  the  smallest.  This  highlights  the  importance  of  perfonning  confirmation 
experiments. 


Table  9.  RPD  Results  for  the  Synthetic  Case  Study 


Robust 

Mean 

Variance 

MSE 

Model 

Point 

Xi 

Goal  — » 

8 

Min 

Min 

Truth 

0.24 

Actual 

7.96 

18.05 

18.05 

Realized 

7.90 

18.63 

18.64 

Std 

0.59 

Predicted 

8.00 

8.62 

8.62 

Realized 

10.82 

19.75 

27.70 

NN 

0.36 

Predicted 

8.00 

24.79 

24.79 

Realized 

9.14 

18.65 

19.95 

CNN 

0.36 

Predicted 

8.00 

24.95 

24.95 

Realized 

9.14 

18.65 

19.95 

CCN 

0.33 

Predicted 

7.81 

19.25 

19.29 

Realized 

8.82 

18.63 

19.30 

KR 

0.30 

Predicted 

8.08 

14.57 

14.58 

Realized 

8.54 

18.62 

18.91 

RBF 

0.27 

Predicted 

7.86 

18.03 

18.05 

Realized 

8.10 

18.62 

18.63 

3.3.2  A  Resistor-Inductor  (RL)  Circuit  Simulation 

An  RPD  study  was  performed  for  a  simulation  of  an  RL  electrical  circuit 
described  by  Kenett  and  Zacks  [82].  The  response  of  interest  is  the  output  current  (in 
amperes)  of  the  circuit  defined  by 
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(55) 


7  =  f/7«2+(2j-/I)1. 

The  four  factors  that  influence  the  output  current  are  listed  in  Table  10.  The  factors  R 
and  L  are  controllable  whereas  the  factors  V  and / are  assumed  to  be  normally  distributed 
random  variables  with  known  means  and  standard  deviations.  The  goal  of  this  RPD 
study  was  to  identify  the  robust  setting  of  R  and  L  that  yields  a  mean  output  current  of  10 
amperes  with  minimum  variation.  A  simulation  of  Equation  (55)  was  created  in 
MATLAB  to  approximate  the  true  mean  and  variance  of  the  system.  To  generate  these 
approximations,  Monte  Carlo  simulations  were  performed  5,000  times  at  900  uniformly- 
spaced  control  design  points.  This  provided  a  reference  point  to  compare  each  modeling 
effort. 


Table  10.  Factors  for  the  RL  Electrical  Circuit  Simulation 


Factor 

Description  (units) 

Min 

Max 

Mean 

Standard 

Deviation 

R 

Resistance  (Q) 

0.05 

9.5 

L 

Self-inductance  (H) 

0.01 

0.03 

V 

Input  voltage  (V) 

100 

3 

f 

Input  frequency  (Hz) 

55 

5/3 

In  this  section,  the  MDM  approach  is  demonstrated  against  two  popular  RPD 
strategies.  Section  3.3.2. 1  compares  the  combined-array  MDM  approach  using  KR  and 
RBF  models  to  the  combined-array  RSM  approach  that  employs  the  NN  model  in 
Equation  (4).  The  NN  model  was  chosen  as  opposed  to  the  Sid  model  in  Equation  (1)  due 
to  the  high  degree  of  non-linearity  in  the  simulation.  Section  3. 3.2. 2  compares  the 
combined-array  RSM  and  MDM  approaches  to  the  combined-array  stochastic  emulator 
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strategy  to  show  that  equivalent  results  can  be  achieved  via  the  former  approach  at  a 
greatly  reduced  computational  cost  and  without  the  need  for  secondary  modeling  efforts. 

S3. 2.1  Combined-Array  MDM  Approach  vs.  Combined-Array  RSM  Approach 

The  NN  model  was  built  from  a  25-run  face-centered  cube  design  which  can  be 
used  within  a  cuboidal  region  to  estimate  the  quadratic  effects  necessary  for  the  response 
model  in  Equation  (4)  [9],  Space-fdling  designs,  such  as  LHS  designs,  are  commonly 
used  for  developing  Kriging  and  RBFNN  models  of  simulations  [83],  Hence,  to  ensure  a 
fair  competition  between  each  model,  the  KR  and  RBF  models  were  each  produced  from 
the  same  25-run  Latin  Hypercube  Sampling  (LHS)  design.  Inspection  of  the  RMSEs  of 
each  model  in  Table  1 1  reveals  that  the  RBF  model  provides  a  slightly  better 
representation  of  the  simulation  than  the  NN  and  KR  models. 

Since  there  are  only  two  control  factors  in  this  RPD  problem,  the  resulting  mean, 
variance,  and  MSE  models  can  be  compared  visually.  These  models  are  shown  in  Figure 
7.  Approximations  of  the  circuit  simulation’s  true  mean,  variance,  and  MSE  are  also 
shown  for  comparison.  The  dots  represent  the  location  of  the  associated  model’s  robust 
point.  Figure  7  shows  that  the  individual  mean  models  provide  the  general  trend  of  the 
simulation’s  true  mean  response.  The  variance  models,  on  the  other  hand,  are  not 
indicative  of  the  simulation’s  true  variance.  The  NN  variance  model  over-estimates  the 
true  variance  throughout  the  entire  design  region.  The  KR  and  RBF  variance  models  are 
better  estimates,  though  they  still  do  not  closely  model  the  non-linearity  of  the  variance 
surface. 
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In  terms  of  robust  points,  NN  and  RBF  provide  similar  points  whereas  KR  most 
closely  identifies  the  location  of  the  true  robust  point.  The  experiment  is  summarized  in 
Table  1 1 .  Each  model’s  robust  point  was  simulated  5,000  times  using  common  random 
numbers  for  V  and/ in  order  to  assess  the  quality  of  the  predicted  values.  All  three 
models  under-estimate  the  simulation’s  mean  at  their  robust  point,  though  the  RBF 
provides  the  closer  prediction.  The  KR  model  offers  the  best  prediction  of  the 
simulation’s  true  variance. 


Figure  7.  Circuit  Simulation  Mean,  Variance,  and  MSE  Models  Using  the  25-run 

Designs 
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Table  11.  Circuit  Simulation  RPD  Results  Using  the  25-run  Designs 


Model 

RMSE 
(%  Range) 

Robust  Point 

R  L 

Mean 

Variance 

MSE 

NN 

1.13 

6.48 

0.019 

Predicted 

9.98 

2.36 

2.36 

(5.5%) 

Realized 

10.73 

0.13 

0.66 

KR 

0.76 

9.10 

0.010 

Predicted 

9.99 

0.12 

0.12 

(4.4%) 

Realized 

10.28 

0.10 

0.18 

RBF 

0.68 

7.15 

0.020 

Predicted 

9.99 

0.35 

0.35 

(3.9%) 

Realized 

10.15 

0.12 

0.14 

Truth 

9.17 

0.011 

10.02 

0.09 

0.09 

At  this  point  in  the  analysis,  it  can  be  concluded  that  after  only  25  experimental 
runs,  the  RSM  and  MDM  approaches  have  provided  adequate,  but  not  highly  accurate, 
representations  of  the  simulation’s  true  mean  and  variance.  However,  since  computer 
simulations  allow  experimenters  to  explore  many  more  factor  levels  and  combinations  of 
factor  levels  than  typically  allowed  in  physical  experiments,  the  number  of  experimental 
runs  in  the  DOE  was  increased  from  25  to  8 1 .  The  NN,  KR ,  and  RBF  models  were  then 
build  from  the  same  81 -run  LHS  design.  Their  corresponding  mean,  variance,  and  MSE 
models  are  depicted  in  Figure  8.  This  new  experiment  is  summarized  in  Table  12.  By 
comparing  Figure  7  and  Figure  8,  improvement  in  the  models’  mean  and  variance 
representations  can  be  observed.  The  KR  and  RBF  mean  models  are  capturing  the 
curvature  of  the  simulation  while  the  NN  mean  model  still  only  represents  its  general 
trend.  There  is  also  significant  improvement  in  the  KR  and  RBF  variance  models.  These 
improvements  can  also  be  seen  in  Table  12  as  the  KR  and  RBF  models’  mean  and 
variance  predictions  at  the  determined  robust  points  are  improved. 
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Figure  8.  Circuit  Simulation  Mean,  Variance,  and  MSE  Models  Using  the  81-run 

LHS  Design 


Table  12.  Circuit  Simulation  RPD  Results  Using  the  81-run  LHS  Design 


Model 

RMSE 
(%  Range) 

Robust  Point 

R  L 

Mean 

Variance 

MSE 

NN 

0.61 

6.60 

0.021 

Predicted 

9.97 

0.52 

0.52 

(2.8%) 

Realized 

10.29 

0.12 

0.20 

KR 

0.37 

2.36 

0.028 

Predicted 

9.98 

0.14 

0.14 

(1.7%) 

Realized 

10.00 

0.17 

0.17 

RBF 

0.34 

8.66 

0.013 

Predicted 

9.98 

0.11 

0.11 

(1.6%) 

Realized 

10.21 

0.10 

0.14 

Truth 

9.17 

0.011 

10.02 

0.09 

0.09 
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Next,  the  number  of  experimental  design  runs  was  increased  further  to  256.  The 
NN,  KR,  and  RBF  models  were  now  built  from  the  same  256-run  LHS  design.  Figure  9 
and  Table  13  show  additional  evidence  that  the  MDM  approach  using  either  the  KR  or 
RBF  models  is  superior  to  the  RSM  approach  using  the  NN  model.  The  mean  and 
variance  surfaces  provide  close  approximations  of  the  simulation’s  true  mean  and 
variance.  The  predicted  means  and  variances  for  each  model  at  their  robust  points  are 
also  nearly  identical  to  their  realized  values.  The  NN  model  never  truly  captures  the  non¬ 
linearity  of  either  the  mean  or  variance  of  the  simulation. 


Figure  9.  Circuit  Simulation  Mean,  Variance,  and  MSE  Models  Using  the  256-run 

LHS  Design 
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Table  13.  Circuit  Simulation  RPD  Results  Using  the  256-run  LHS  Design 


Model 

RMSE 
(%  Range) 

Robust  Point 

R  L 

Mean 

Variance 

MSE 

NN 

0.49 

6.52 

0.021 

Predicted 

9.98 

0.43 

0.43 

(2.8%) 

Realized 

10.24 

0.12 

0.18 

KR 

0.05 

8.46 

0.016 

Predicted 

9.99 

0.08 

0.08 

(0.3%) 

Realized 

9.99 

0.10 

0.10 

RBF 

0.08 

7.34 

0.020 

Predicted 

9.99 

0.11 

0.11 

(0.4%) 

Realized 

9.99 

0.11 

0.11 

Truth 

9.17 

0.011 

10.02 

0.09 

0.09 

3. 3. 2. 2  Combined-Array  RSM/MDM Approaches  vs.  Combined-Array 
Emulator  Approach 

The  combined-array  RSM  and  MDM  approaches  were  compared  to  the 
combined-array  stochastic  emulator  approach  using  the  NN,  KR,  and  RBF  models 
generated  via  the  256-run  LHS  design  as  the  emulators.  Within  the  emulator  approach, 
the  individual  mean  and  variance  models  were  built  from  a  5"  factorial  design  (labeled 
DOE2  in  Figure  4).  These  models  were  constructed  using  the  same  modeling  approach 
that  created  the  emulator  itself.  That  is,  if  Kriging  was  used  to  build  the  emulator,  then 
Kriging  was  also  used  to  build  the  associated  mean  and  variance  models. 

The  emulator  strategy  requires  a  large  number  of  replications  to  be  performed  at 
each  design  point.  To  determine  an  appropriate  number  of  replications,  the  long-run 
mean  and  variance  of  each  model  was  examined  at  the  center  point  and  the  four  corner 
points  of  the  design  (Figure  10).  Based  on  these  plots,  the  choice  was  made  to  run  500 
replications  at  each  of  the  25  design  points  of  DOE2  since  both  the  mean  and  variance  of 
each  model  tend  to  reach  a  steady-state  at  that  point.  Therefore,  the  emulator  approach 
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required  25  x  500  =  12, 500  evaluations  of  the  emulator  in  order  to  build  the  mean  and 


variance  models.  This  is  in  addition  to  the  original  256  experiments  necessary  to  build 
the  emulator  in  the  first  place. 
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Figure  10.  Long-Run  Analysis  of  NN,  KR,  and  RBF  Circuit  Simulation  Emulators 


A  comparison  of  the  RSM  and  MDM’s  mean,  variance,  and  MSE  models  in 
Figure  9  to  the  emulator’s  corresponding  models  in  Figure  1 1  reveals  very  similar 
response  surfaces.  A  further  examination  of  the  robust  point  summaries  for  the 
RSM/MDM  (Table  13)  and  emulator  (Table  14)  approaches  also  shows  nearly  equivalent 
results.  The  robust  points,  as  well  as  their  predicted/realized  means  and  variances,  for 
each  model  type  are  approximately  equal.  Consequently,  it  can  be  concluded  that  the 
combined-array  RSM  and  MDM  approaches  and  the  stochastic  emulator  strategy  yielded 
similar  results.  However,  the  emulator  strategy  required  the  development  of  two 
experimental  designs,  three  response  modeling  efforts  (one  each  for  the  emulator,  the 
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mean,  and  the  variance),  one  optimization  procedure,  and  12,500  additional  function 
evaluations.  The  RSM/MDM  approach,  on  the  other  hand,  only  required  one 
experimental  design,  one  response  modeling  effort,  and  one  optimization  procedure. 


Figure  11.  Emulator  Approach  Mean,  Variance,  and  MSE  Models  for  the  Circuit 

Simulation 
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Table  14.  Emulator  Approach  Results  for  the  Circuit  Simulation 


Rnhust  Point 

Variance 

MSE 

Model 

R 

L 

Mean 

NN 

6.45 

0.021 

Predicted 

9.98 

0.39 

0.39 

Realized 

10.22 

0.12 

0.17 

KR 

8.32 

0.016 

Predicted 

9.99 

0.08 

0.08 

Realized 

10.06 

0.10 

0.10 

RBF 

7.33 

0.020 

Predicted 

9.99 

0.10 

0.10 

Realized 

9.97 

0.11 

0.11 

Truth 

9.17 

0.011 

10.02 

0.09 

0.09 

3.4  Summary 

Chapter  III  has  extended  the  combined-array  RSM  approach  that  relies  upon  low- 
order  polynomial  models.  It  was  demonstrated  that  improved  models  of  a  computer 
simulation’s  mean  and  variance  can  be  attained  through  the  MDM  approach  that  employs 
non-linear  response  modeling  techniques  such  as  Kriging  or  RBFNN  models.  It  was 
further  shown  that  the  combined-array  MDM  approach  generates  results  that  are 
approximately  equivalent  to  the  stochastic  emulator  approach.  However,  these  results 
can  be  achieved  at  a  greatly  reduced  computational  cost  by  utilizing  the  MDM  approach. 

Throughout  this  chapter,  the  concern  was  more  with  examining  the  individual 
mean  and  variance  models  and  less  with  the  actual  location  of  the  robust  points.  In  fact, 
each  of  the  12  robust  points  that  were  identified  for  the  circuit  simulation  are  very  good 
solutions.  This  is  illustrated  in  Figure  12.  The  points  represented  with  are  design 
space  locations  in  which  the  true  MSE  of  the  simulation  is  less  than  or  equal  to  0.25. 
Actually,  these  points  represent  a  region  of  points  and  not  just  individual  points.  The 
minimum  MSE,  denoted  with  “X”,  is  0.09.  By  comparison,  the  maximum  MSE  in  the 
region  is  361.18.  The  12  individual  robust  points  are  denoted  with  “O”.  It  is  shown  that 
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they  each  fall  in  or  near  a  “robust  region”  which  implies  that  they  are,  in  their  own  right, 
robust  solutions  for  this  RPD  problem.  The  actual  interest,  in  this  case,  lies  with  the 
predicted  means  and  variances  at  these  robust  points.  It  was  shown  that  the  MDM 
approach  was  superior  to  the  RSM  approach  in  providing  improved  predictions  of  the 
system. 

0.03 

0.025 

-J  0.02 

0.015 

0.01 

Figure  12.  Illustration  of  the  Circuit  Simulation’s  “Robust  Region” 
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IV.  A  Nested  Desirability  Function-Based  Approach  to  Multi-Response  RPD 


4.1  Introduction 

This  chapter  considers  the  situation  in  which  an  experimenter  seeks  to  examine 
the  influence  a  set  of  independent  variables  has  on  several  system  responses 
simultaneously.  For  example,  a  situation  may  require  finding  the  optimal  set  of 
conditions  that  reduce  cost  while  also  increasing  yield.  Unfortunately,  the  solution  that 
minimizes  cost  is  most  likely  not  the  same  as  the  solution  that  maximizes  yield.  In  fact, 
the  two  solutions  may  diametrically  oppose  each  other.  In  a  case  such  as  this,  trade-offs 
between  the  responses  must  be  explored  to  find  a  collaborative  solution. 

Chapter  IV  focuses  on  the  second  objective  of  this  dissertation:  to  develop  an 
approach  for  multi-response  RPD  problems  that  provides  a  collaborative  solution  that  is 
balanced  across  the  means  and  variances  of  each  response.  Existing  techniques  seek  to 
find  an  optimal  balance  across  the  set  of  responses.  However,  in  some  cases  the  mean  or 
variance  of  one  response  influences  the  solution  in  such  a  way  that  the  means  and 
variances  of  the  remaining  responses  are  insignificant  to  the  overall  RPD  problem.  In 
these  situations,  it  is  difficult  to  attain  a  solution  that  is  balanced  across  all  responses. 

This  chapter  presents  a  technique  based  on  well-known  desirability  functions  that 
places  the  means  and  variances  of  all  responses  on  a  level  playing  field.  The  proposal 
also  allows  decision  makers  to  integrate  their  personal  preferences  for  the  individual 
response  means  and  variances.  For  example,  they  may  be  willing  to  accept  sacrificing 
being  slightly  “off-target”  in  one  response,  but  be  reluctant  to  allow  the  variation  of 
another  response  to  get  too  large. 
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The  rest  of  this  chapter  continues  in  the  following  manner.  Section  4.2  describes 
the  proposed  desirability  function-based  methodology  for  solving  multi-response  RPD 
problems.  Section  4.3  utilizes  two  case  studies  to  compare  the  proposed  approach  to  a 
popular  MSE-based  approach.  Finally,  Section  4.4  summarizes  this  chapter. 

4.2  The  Nested  Desirability  Function  Approach 

The  proposed  approach  to  solving  the  multi-response  RPD  problem  employs  the 
desirability  functions  established  by  Derringer  and  Suich  [14].  However,  instead  of  using 
them  to  transform  each  response  into  a  corresponding  desirability  level,  they  were  used  to 
transform  the  individual  response  means  and  variances  into  desirability  levels.  This 
approach  allows  a  decision  maker  to  state  their  preferences  regarding  what  is,  and  is  not, 
acceptable  for  each  response’s  mean  and  variance.  Park  and  Park’s  [17]  weighted 
desirability  function  approach  for  optimization  problems  was  also  implemented  to  apply 
distinct  degrees  of  importance  to  the  different  responses,  response  means,  and  response 
variances. 

Let  j>;  denote  an  estimated  relationship  between  response  variable  j  and  the  vector 
of  independent  variables  x.  Also,  suppose  that  juj  and  cr  are  the  respective  mean  and 
variance  of  j>. .  The  desirability  of  response  j  is  defined  as 

<56> 

where  0<djM<\  and  0  <d  ,  <  1  are  the  respective  desirability  transfonnations  for  //,  and 
av  .  Also,  Wj  M>  0  and  w  .  >  0  are  the  preferred  weights  for  juj  and  ov  ,  respectively. 

Since  there  are  two  desirabilities — one  for  the  mean  and  one  for  the  variance — the 
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relationship  w.  /(  +  w.  =  2  must  hold.  D; ,  which  is  contained  in  the  interval  [0,1],  will 

be  0  if  either  d  j  or  d  are  0.  This  implies  that  response  j  is  unacceptable  if  either  its 

mean  or  variance  is  unacceptable.  Finally,  the  overall  desirability  of  the  combined  mean 
and  variance  levels  across  all  J  responses  is 


D  = 


(57) 


where  W t  >  0  is  the  preferred  weight  for  response  j  and  ^  IF.  -  J  .  Again,  if  Dj  =  0  for 

M 

any  j,  then  D  =  0  and  the  whole  product  is  unacceptable.  The  robust  point  is  the  x  that 
maximizes  D.  Since  D  is  composed  of  the  desirabilities  Dl,D^,...,DJ  and  each  D/  is  itself 

composed  of  the  desirability  functions  d u  and  d .  ^  ,  the  proposed  procedure  will  be 
referred  to  as  the  Nested  Desirability  ( ND )  approach. 


4.2.1  Desirability  Transformations  for  Juj 

Taguchi  [4]  specified  three  experimental  cases  for  managing  a  system’s  mean 
response  in  a  RPD  scenario:  smaller-the-better,  larger-the-better,  and  target-is-best.  Let 
Lj  and  Uj  be  the  respective  minimum  and  maximum  values  of  / / , .  If  response  j  is  a 

smaller-the-better  case,  then  the  desirability  transformation  for  ju  is 


dJ.r=< 


r  U.  ~jU  ^ 

J4‘ 

V  Uj,M  ~  j 


L  <  Li  <U 

J.M  O  J-U 

otherwise 


(58) 


If  response  j  is  a  larger-the-better  case,  then  the  desirability  transformation  for  jj .  is 
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dJ.M=< 


(  A  J  V 

U  -L 

“j  J'M 

KUJ*-Li*J 


L  <  u.<U ■ 

J,M  “ J  J-M 

otherwise 


(59) 


Finally,  if  response  j  is  a  target-is-best  case  with  a  desired  target  of  r  where 


L  t  ,  then  the  desirability  transfonnation  for  ju .  is 


r  ~  T  \ 

thzhi 

\xrLi.n 

f  TT  ->  A 

Uj*-» 


r/<A <L,, 


(60) 


I  0  otherwise 

The  exponents  r  and  5  in  Equations  (58)-(60)  are  shape  parameters  for  the  desirability 
function.  A  value  for  r  or  5  that  is  greater  than  1  puts  greater  importance  on  the  mean 
response  values  being  closer  to  the  respective  minimum,  maximum,  or  target  value.  On 
the  other  hand,  a  value  for  r  or  s  that  is  between  0  and  1  implies  that  the  desirability  value 
is  large  even  if  the  mean  response  is  far  from  its  goal  value. 


4.2.2  Desirability  Transformations  fora2 

Whereas  there  are  three  distinct  cases  for  managing  jjj ,  there  is  only  one  such 
case  fora2  since  a  minimum  response  variance  is  always  desired.  Therefore,  letZ  ^  and 
U  n  be  the  respective  minimum  and  maximum  values  of  a2 .  Then,  the  desirability 
transformation  fora2  is 
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(61) 


d  2  =  < 


U.  2  -of  A 

J.O  J 

U  2-L  4 

V  y'.CT  j,<r  J 


L  2  <a2<U  2 

y.o"  7  y.o-- 


otherwise 


The  value  r  is  defined  as  it  was  in  Section  4.2. 1 . 


4.3  Application  and  Results 

This  section  demonstrates  Koksoy  and  Yalcinoz’s  (AT)  MSE  procedure  [30]  and 
the  Nested  Desirability  (ND)  procedure  on  two  case  studies:  a  synthetic  case  consisting  of 
two  known  functions  and  a  textbook  example  for  a  physical  experiment. 

4.3.1  A  Synthetic  Case  Study 

The  response  functions  in  Equations  (61)  and  (62)  were  used  to  demonstrate  the 
ND  procedure  alongside  the  KY  procedure  in  Equation  (27).  The  responses,  which  were 
considered  to  be  equally  important,  were  influenced  by  the  control  factor  x  e  [-1,1]  and 

the  noise  factor  z  ~  N(jUz  =  0,  a2  =  l) . 

yl  =  -4x2  -3x  +  4z2  +3z  +  4xz  +  15  (62) 

y2  =  -jx-jz2 -jz  +  jxz  +  6  (63) 

Given  the  known  distributional  parameters  for  z,  the  means  and  variances  of  Equations 

(62)  and  (63)  can  be  detennined  analytically.  These  functions  are  shown  in  Equations 

(63) -(66). 


//,  =  -4x2  -3x  +  19 

(64) 

of  =  (4x  +  3)2  +32 

(65) 
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(66) 


Pi  ^  x  +  2 

^2=(f^-})2+l  (67) 

In  this  scenario,  y]  was  a  larger-the-better  response  and  v2  was  a  smaller-the-better 
response. 

4.3. 1.1  The  Koksoy  and  Yalcinoz  (KY)  Procedure 
Since  the  responses  were  deemed  equally  important,  the  weights  in  the  KY 
procedure  were  W1  =  W2  =  \ ,  Thus,  the  jointly  robust  point  xK}  was  found  by  minimizing 

KY  =  jMSEl  +jMSE2  where  MSEX  =  ay  -//,2  and  MSE2  =  o\  +  ji\ .  Also,  in  order  to 
examine  any  trade-offs  made  among  the  responses  in  detennining  a  jointly  robust 
operating  point,  xA>  was  compared  to  the  results  of  the  two  single  response  RPD 
problems  suggested  by  Equations  (62)  and  (63).  That  is,  the  robust  point  for^  was 

identified  without  any  consideration  ofy2  by  minimizing  MSE]  to  yield  the  solution xfy  . 
Similarly,  MSE2  was  minimized  to  find  the  robust  point  forjy  only.  This  solution  is 
identified  as  xf Y . 

Figure  13  illustrates  the  results  of  solving  the  three  RPD  problems.  The  locations 
of  the  single  response  robust  solutions  are  shown  in  Figure  13a  and  Figure  13b.  Clearly, 
these  individual  solutions  occur  at  opposite  ends  of  the  range  of  x  suggesting  that  some 
compromise  will  need  to  be  made  to  identify  an  operating  point  that  is  mutually  robust 
for  y,  and  y2 .  Finally,  the  jointly  robust  pointxA>  is  shown  in  Figure  13c.  It  can  be 

observed  that,  even  though  the  responses  were  given  equal  weights,  xA>  is  nearly 
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identical  to  xf Y .  It  appears  as  if  there  was  no  compromise  amongst  the  responses.  This 


can  be  further  examined  by  decomposing  KY. 


MSE1  MSE2 


xxx 
(a)  (b)  (c) 

Figure  13.  Robust  Point  Locations  for  the  Synthetic  Case  Using  the  KY  Procedure 

Figure  14  plots  KY  with  its  decomposed  functions,  namely \MSE{  and \MSE2 ,  on 
the  same  graph.  It  is  now  clear  thatxAi  is  indistinguishable  from  due  to  the  fact  that 
\MSE2  is  relatively  unchanging  when  compared  to  \  MSEX  and  has  little  to  no  effect  on 
KY.  The  range  of \MSEl  is  142.8  while  the  range  of  \MSE2  is  only  8.95.  Based  on  this 
simple  decomposition,  it  can  be  concluded  that  the  jointly  robust  solution  is  greatly 
influenced  by  y1  whereas  y2  has  a  negligible  effect  on  the  joint  solution. 
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Figure  14.  Decomposition  of  KY  for  the  Synthetic  Case  Study 


Each  response’s  MSE  function  can  be  further  examined  by  decomposing  them 
into  their  variance  and  squared  bias  terms.  In  the  case  ofMSEj ,  its  variance  and  negative 

squared  bias  term  (-juf)  will  actually  be  investigated.  Again,  it  is  shown  in  Figure  15a 

and  Figure  15b  that  the  MSE  functions  are  driven  by  a  single  term — in  each  case,  the 
squared  bias.  The  variance  terms  have  minimal  influence  on  the  MSE  functions  and,  as  a 
result,  minimal  influence  on  determining  xf 5  ,  xf1  ,  or  xAi  .  It  can  be  concluded  that,  for 
this  example,  the  jointly  robust  solution  is  influenced  significantly  by  //, .  Conversely,  ju2 
,  of  ,  and  of  have  limited,  if  any,  impact  on  the  overall  RPD  solution  derived  via  the  KY 
procedure. 
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Figure  15.  Decomposition  of  Individual  MSE  Functions  for  the  Synthetic  Case 

Study 


43.1.2  The  Nested  Desirability  (ND)  Procedure 

Linear  (r  =  1)  desirability  functions  were  utilized  for  the  ND  procedure.  Again, 
since  the  responses  hold  equal  importance,  the  weights  in  the  ND  formulation  of  the 
problem  were  Wx  =  W2  =  1 .  Now,  the  KY  procedure  does  not  consider  weighting  the 
contributions  of  the  individual  means  and  variances;  here,  the  ND  procedure  weighted 
them  equally.  Specifically,  wj  fJ  =  w.  ^  - 1  for  j  =  1,2  .  The  jointly  robust  point  x™ 

maximizes  D  =  ((D1)1  x(D2)1)/^  where  D{  =  x(dla2 2  and 


D2  =  j'  x  (cL,  ^  j'  j  .  Even  though  the  weights  are  equal  to  1 ,  they  were  included 

here  as  the  exponents  for  completeness.  Again,  the  results  of  the  two  single  response 
RPD  problems  are  presented.  That  is,  the  robust  point  for  y! ,  denoted  x,™  ,  was  found  by 
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maximizing  D\.  Similarly,  FT  was  maximized  to  find  the  robust  point  fory,  only.  This 
solution  is  identified  as  . 

Figure  16  illustrates  the  results  of  solving  the  three  RPD  problems  using  the  ND 
approach.  The  locations  of  the  single  response  robust  solutions  are  shown  in  Figure  16a 
and  Figure  16b.  Again,  as  the  KY procedure  showed,  these  solutions  are  conflicting. 
Finally,  the  jointly  robust  point  x  v/'  is  shown  in  Figure  16c.  As  opposed  to  the  results  of 
the  KY  procedure,  some  compromise  can  now  be  observed  between  the  two  responses  in 
order  to  generate  a  jointly  robust  point.  Table  15  and  Table  16  can  be  examined  further 
in  Section  4.3. 1.3  to  see  where  these  compromises  were  made. 


(a) 


(b) 


(c) 


Figure  16.  Robust  Point  Locations  for  the  Synthetic  Case  Using  the  ND  Procedure 


43.1.3  Comparison  of  the  KY  Procedure  and  the  ND  Procedure 
The  results  of  the  KY  procedure  and  the  ND  procedure  can  be  observed  through 
two  different  points-of-view  (POV).  Table  15  summarizes  the  means,  variances,  and 
MSEs  of  the  three  robust  points  for  the  two  RPD  strategies  by  taking  a  “MSE  POV.”  The 
italicized  values  are  the  corresponding  results  of  the  secondary  response  when  each  single 
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response  RPD  problem  is  solved.  The  bold  values  are  the  results  of  solving  the  multi¬ 
response  RPD  problems.  As  noted  in  Figure  13a  and  Figure  13c,  x,A>  and.rAI  are  nearly 
equal.  Therefore,  when  an  equally  important  y2  is  added  into  the  problem  formulation, 
there  exists  a  slight  trade-off  by  increasing  of  for  marginal  improvements  in  //2  and  crl . 

For  all  intents  and  purposes,  there  is  practically  zero  compromise  between  the  responses 
in  the  joint  RPD  case  when  using  the  KY  procedure.  This  is  not,  however,  the  situation 
when  using  the  ND  procedure  where  noticeable  trade-offs  occur  amongst  the  means  and 
variances  of  the  two  responses.  Now,  when  y2  is  considered  as  important  as  yl ,  larger 
degradations  in  //,  and  of  are  incurred  in  exchange  for  substantial  improvements  in  ju2  and 
of  .  Based  on  the  decision  maker’s  original  preferences,  this  is  a  solution  they  are 

willing  to  accept.  However,  by  observing  the  results  from  an  MSE  POV,  xND  represents 
a  12.51%  decrease  in  overall  “ KY  value”  when  compared  to  xKY .  This  is  understandable 
since  xKi  is  the  best  solution  in  terms  of  MSE. 
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Table  15.  Synthetic  Multi-Response  Case  Study  from  a  MSE  Point-of-View 


Robust 

Point 

RPD  Goal  — * 

Mi 

Max 

Min 

MSE l 

Min 

Mi 

Min 

Min 

mse2 

Min 

KY 

Min 

4Y  = 

-0.41 

19.56 

33.85 

-348.74 

5.83 

0.63 

34.62 

-157.06 

xf  = 

=  1.00 

12.00 

81.00 

-63.00 

4.70 

0.54 

22.63 

-20.19 

.vA1  = 

:  -0.38 

19.56 

34.17 

-348.42 

5.81 

0.62 

34.38 

-157.02 

_ 

X\  ~ 

-0.52 

19.48 

32.86 

-346.61 

5.91 

0.67 

35.60 

-155.51 

_ 

xi  - 

=  1.00 

12.01 

80.93 

-63.31 

4.70 

0.54 

22.63 

-20.34 

xND  - 

=  0.10 

18.66 

43.53 

-304.67 

5.42 

0.53 

29.91 

-137.38 

%  Improvement 
from  xK)  tox* 

-4.60% 

-27.39% 

-12.56% 

6.71% 

14.52% 

13.02% 

-12.51% 

Additionally,  the  results  can  also  be  examined  from  a  “desirability  POV.”  Table 
16  summarizes  the  desirabilities  for  the  means,  variances,  and  responses  of  the  three 
robust  points  for  the  two  RPD  strategies.  The  trade-offs  in  the  desirability  POV  are 
similar  to  the  trade-offs  in  the  MSE  POV.  However,  when  looking  at  the  same  results 
through  a  desirability  lens,  utilizing  xND  as  the  jointly  robust  point  instead  ofxA>  results  in 
12.00%  and  20.83%  decreases  in  desirability  for  //,  and  of ,  respectively.  In  exchange  for 
this,  however,  the  desirability  levels  increase  significantly — 77.42%  and  43.08% — for //2 
and  of  .  This  ultimately  results  in  a  16.67%  increase  in  the  overall  system  desirability  at 
xND  .  Again,  this  is  expected  since  xM>  is  the  best  solution  in  terms  of  desirability. 
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Table  16.  Synthetic  Multi-Response  Case  Study  from  a  Desirability  Point-of-View 


Robust 

Point 

D\ 

4, 

Di 

D 

xfV=- 

-0.41 

1.00 

0.96 

0.98 

0.30 

0.63 

0.43 

0.65 

4¥  = 

1.00 

0.00 

0.00 

0.00 

1.00 

0.89 

0.94 

0.00 

xKY  =  ■ 

-0.38 

1.00 

0.96 

0.98 

0.31 

0.65 

0.45 

0.66 

_ 

X1  — 

-0.52 

0.99 

0.98 

0.99 

0.24 

0.54 

0.36 

0.60 

_ 

X2  ~ 

1.00 

0.00 

0.00 

0.00 

1.00 

0.89 

0.94 

0.04 

xND  = 

0.10 

0.88 

0.76 

0.82 

0.55 

0.93 

0.71 

0.77 

%  Improvement 
fromxA1  to  xNn 

-12.00% 

-20.83% 

-16.33% 

77.42% 

43.08% 

57.78% 

16.67% 

4.3.2  The  Force  Transducer  Experiment 

Details  for  the  following  example  can  be  found  in  Romano  et  al.  [22].  Several 
authors  have  also  used  this  example  as  a  case  study  for  their  multi-response  RPD 
methods  [30-32].  In  short,  the  problem  consists  of  two  response  variables  (y,  and  v2) , 

three  control  variables  (x15x2  and  x3) ,  and  two  noise  variables  (z1  and  z2)  .  The 
experimental  results  are  displayed  in  Table  17.  The  noise  factors  were  assumed  to  be 
independent  with  zero  mean  and  variances  of  and  of  .  Thus,  the  ±1  experimental  levels 
for  z j  were  set  to  ±cr  . .  Per  Romano  et  al.  [22],  the  fitted  response  surface  functions  for  Vj 
and  v2  are 

j>j(x,z)  =  1.38-0. 361xj  -0.155x2  +  0.0771x3  -0.148x,x2  +0.0218X[X3 
+0.0130x2x3  +0.0481xf  -0.0588zj  -0.01 16z2  -1-O.OIOOxjZj 

and 
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j>2(x,z)  =  1.64-0. 592xj  +0.438x2  -0.0950x3  +0.301xjX2  -0.143xjX3 
+0.20  lx2  -0.0844xjX2x3  +  0.0794x1z1 


(69) 


The  authors  reported  that  no  lack  of  fit  was  detected  and  the  MSEs  of  the  models  were 
a 2,  =  0.0003253  and  <t22  =  0.024  . 


Table  17.  Experimental  Results  for  the  Force  Transducer  Experiment 


Run 

*2 

*3 

Zi 

^2 

Ti 

y2 

1 

-1 

-1 

-1 

-1 

1 

1.81 

1.10 

2 

-1 

-1 

-1 

1 

-1 

1.69 

l.n 

3 

-1 

-1 

1 

-1 

-1 

1.90 

1.07 

4 

-1 

-1 

1 

1 

1 

1.78 

1.07 

5 

-1 

1 

-1 

-1 

-1 

1.80 

1.47 

6 

-1 

1 

-1 

1 

1 

1.63 

1.18 

7 

-1 

1 

1 

-1 

1 

1.92 

1.41 

8 

-1 

1 

1 

1 

-1 

1.78 

1.58 

9 

1 

-1 

-1 

-1 

-1 

1.36 

1.57 

10 

1 

-1 

-1 

1 

1 

1.22 

2.03 

11 

1 

-1 

1 

-1 

1 

1.48 

1.38 

12 

1 

-1 

1 

1 

-1 

1.44 

1.68 

13 

1 

1 

-1 

-1 

1 

0.693 

3.37 

14 

1 

1 

-1 

1 

-1 

0.616 

3.75 

15 

1 

1 

1 

-1 

-1 

0.950 

2.81 

16 

1 

1 

1 

1 

1 

0.817 

2.83 

17 

-1 

0 

0 

0 

0 

1.79 

1.24 

18 

1 

0 

0 

0 

0 

1.03 

2.46 

19 

0 

-1 

0 

0 

0 

1.53 

1.23 

20 

0 

1 

0 

0 

0 

1.22 

1.73 

21 

0 

0 

-1 

0 

0 

1.30 

1.63 

22 

0 

0 

1 

0 

0 

1.44 

1.67 

23 

0 

0 

0 

0 

0 

1.38 

1.73 

24 

0 

0 

0 

0 

0 

1.39 

1.74 

25 

0 

0 

0 

0 

0 

1.40 

1.74 

The  goal  of  this  RPD  study  was  to  find  the  robust  settings  ofx1?x2  and  x3  that 


minimized yx  andy2 .  Similar  to  the  synthetic  case  study  in  Section  4.3.1,  this  problem 
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was  solved  using  the  KY procedure  and  the  ND  procedure.  The  estimated  mean  models 


are 

ju{  =  1.38-0. 361xj  -0.155x2  +  0.0771x3  -0.148XjX2  +  0.0218XjX3 
+0.0130x2x3  +  0.048  lxf 

and 

ju2  =  1.64-0. 592xj  +  0.438x2  -0.0950x3  +0.301XjX2  -0.143x!X3 
+0.201xj2  -0.0844X[X2x3 

The  estimated  variance  models  are 

<7?  =  (-0.0588  +  O.OlXj  )2  +  (-0.01 16)2  +  <j2el 


and 


<r2  =  (0.0794xj  )2  +a22. 


(70) 


(71) 


(72) 


(73) 


Table  18  summarizes  the  means,  variances,  and  MSEs  of  the  three  robust  points 
for  the  two  RPD  strategies  from  a  “MSE  POV.”  In  this  case,  as  opposed  to  the  synthetic 
case  study,  both  the  KY  procedure  and  the  ND  procedure  generated  solutions  in  which 
compromises  were  made  amongst  the  means  and  variances.  However,  xND  represents  a 
33.29%  decrease  in  overall  “ KY  value”  when  compared  to  xA  >  .  Similar  conclusions 
regarding  trade-offs  between  the  response  means  and  variances  can  be  made  by 
examining  Table  19.  Now,  from  a  “desirability  POV,”  utilizingxv/)  as  the  jointly  robust 
point  instead  ofxA1  yields  a  15.25%  increase  in  the  overall  system  desirability. 
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Table  18.  Force  Transducer  Experiment  Results  from  a  MSE  Point-of-View 


Robust  Point 

Mx 

MSEs 

Mi 

07 

mse2 

KY 

RPD  Goal  — ► 

Min 

Min 

Min 

Min 

Min 

Min 

Min 

xf7  =(1.00,1.00,-1.00) 

0.65 

0.0028 

0.43 

3.49 

0.0302 

12.21 

6.32 

xf  =(-0.57,-1.00,1.00) 

1.72 

0.0046 

2.96 

1.04 

0.0260 

1.11 

2.04 

xKY  =(0.07,-1.00,1.00) 

1.59 

0.0038 

2.53 

1.12 

0.0240 

1.28 

1.91 

x™  =(1.00,1.00,-1.00) 

0.65 

0.0028 

0.43 

3.49 

0.0302 

12.21 

6.32 

xf5  =(-0.04,-1.00,1.00) 

1.61 

0.0040 

2.60 

1.10 

0.0240 

1.23 

1.92 

xND  =(0.32,-0.20,-1.00) 

1.23 

0.0036 

1.52 

1.88 

0.0246 

3.56 

2.54 

%  Improvement 

22.64% 

5.26% 

40.29% 

-67.86% 

-2.50% 

-176.79% 

-33.29% 

Table  19.  Force  Transducer  Experiment  Results  from  a  Desirability  Point-of-View 


Robust  Point 

Kn 

Di 

D1 

D 

xf7  =(1.00,1.00,-1.00) 

1.00 

1.00 

1.00 

0.00 

0.00 

0.00 

0.00 

xfY  =(-0.57,-1.00,1.00) 

0.11 

0.24 

0.16 

1.00 

0.68 

0.82 

0.36 

xKY  =(0.07,-1.00,1.00) 

0.22 

0.58 

0.36 

0.97 

1.00 

0.98 

0.59 

x(®  =(1.00,1.00,-1.00) 

1.00 

1.00 

1.00 

0.00 

0.00 

0.00 

0.00 

xf3  =(-0.04,-1.00,1.00) 

0.20 

0.52 

0.32 

0.98 

1.00 

0.99 

0.57 

xND  =(0.32,-0.20,-1.00) 

0.52 

0.70 

0.60 

0.66 

0.90 

0.77 

0.68 

%  Improvement 
from  xKY  to  xND 

136.36% 

20.69% 

66.67% 

-31.96% 

-10.00% 

-21.43% 

15.25% 

4.4  Summary 

MSE-based  strategies  for  solving  multi-response  RPD  problems  are  popular 
methods.  However,  as  the  synthetic  case  study  showed,  sometimes  the  mean  or  variance 
of  one  response  can  render  the  means  and  variances  of  the  remaining  responses 
insignificant  to  the  overall  RPD  problem.  This  can  diminish  the  chance  of  finding  a 
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compromised  solution.  Also,  trade-offs  among  the  means  and  variances  of  the  individual 
responses  are  typically  ignored.  Chapter  IV  presented  an  approach,  based  on  well-known 
desirability  functions,  that  is  beneficial  in  two  ways.  First,  it  places  the  responses,  as  well 
as  their  means  and  variances,  on  equal  footing.  Second,  it  allows  a  decision  maker  to 
declare  their  personal  preferences  for  the  responses’  means  and  variances  from  the  outset. 
The  resultant  operating  point  is  a  system  setting  that,  whether  one  is  looking  at  the 
problem  from  a  MSE  POV  or  a  desirability  POV,  produces  a  mutually  robust  set  of 
responses  the  decision  maker  considers  acceptable. 
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V.  Quality  Measures  for  Comparing  Multiple  RPD  Strategies 


5.1  Introduction 

There  is  a  growing  literature  in  which  multiple  RPD  problem  solving  strategies 
are  contrasted.  Articles  typically  report  a  system’s  predicted  mean  and  variance  found 
using  a  number  of  approaches;  however,  some  neglect  to  demonstrate  the  quality  of  those 
predictions  through  confirmation  experiments  [22,  24,  25,  29-32,  46,  71].  There  are  two 
operative  questions.  Are  the  predictions  good  estimates  of  the  system’s  actual 
perfonnance  at  the  estimated  robust  settings?  Also,  which  method’s  robust  point 
estimate  actually  realizes  the  best  combination  of  mean  and  variance?  These  questions 
become  more  challenging  when  there  are  multiple  responses  of  interest.  To  provide  a 
framework  for  addressing  these  questions,  this  research  first  posits  that  any  RPD  problem 
solving  methodology  can  be  evaluated  on  the  basis  of  the  accuracy  of  its  predictions  and 
its  realized  robustness  relative  to  its  competitors. 

Chapter  V  addresses  the  third  objective  of  this  dissertation  which  is  aimed  at 
generating  measures  for  comparing  different  RPD  problem  solving  strategies.  Such 
measures  can  increase  the  understanding  of  each  strategy  and  allow  the  analyst  to  make  a 
more  knowledgeable  evaluation  of  the  competing  procedures.  The  rest  of  Chapter  V  is 
organized  as  follows.  Section  5.2  introduces  a  methodology  for  meeting  this  research 
objective.  Section  5.3  conducts  a  case  study  using  a  discrete  event  simulation  that 
compares  the  results  of  12  competing  robust  design  strategies.  Finally,  Section  5.4 
summarizes  this  chapter. 
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5.2  Quality  Indices  for  Robust  Parameter  Design 

As  mentioned,  the  RPD  literature  generally  compares  the  system’s  predicted 
mean  and  variance  at  the  robust  setting  detennined  through  one  approach  to  the  system’s 
predicted  mean  and  variance  at  the  robust  setting  found  though  a  competing  method. 
Authors  then  make  a  choice  about  which  method  is  superior  based  upon  the  closeness  of 
these  predictions  to  ideal  mean  and  variance  targets.  For  example,  in  Taguchi’s  target-is- 
best  scenario,  one  would  prefer  a  robust  point  with  a  predicted  mean  close  to  the  target 
value  with  a  small  predicted  variance  over  another  setting  with  a  predicted  mean  further 
from  the  target  having  a  higher  predicted  variance.  Based  on  the  predictions,  it  is 
expected  that  the  fonner  operating  point  is  better  than  the  latter.  It  is  this  comparison  that 
drives  some  authors  to  conclude  the  dominance  of  one  methodology  over  another. 
However,  if  the  methodologies  provide  poor  predictions,  then  their  conclusions  may  be 
inaccurate. 

To  make  a  more  comprehensive  assessment  of  the  competing  strategies,  this 
research  recommends  expanding  the  analysis  to  include  confirmation  experiments  at  each 
robust  setting.  First,  it  must  be  stated  that  this  is  obviously  not  feasible  in  all  cases.  It  is 
highly  unlikely  that  a  capital-producing  manufacturing  line  would  shut  down  in  order  to 
make  some  trial  runs.  It  may  also  be  extremely  costly  or  dangerous  to  make  these  test 
runs  in  other  situations.  However,  these  confirmatory  runs  may  be  easy  and  inexpensive 
to  perfonn  when  simulations  are  used.  By  performing  confirmation  experiments  at  each 
robust  setting,  further  insight  can  be  gained  into  each  RPD  approach’s  accuracy  and 
robustness  qualities.  Accuracy  is  a  quality  measure  of  how  close  the  predicted  mean  and 
variance  of  the  system  are  to  the  realized  mean  and  variance  of  the  system  when  it  is 
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repeatedly  executed  at  a  specific  operating  point.  Robustness,  on  the  other  hand,  assesses 
how  close  the  system’s  realized  mean  and  variance  are  to  their  desired  target  values.  A 
minimum  response  variance  is  always  preferred;  however,  Taguchi’s  overall 
experimental  goal  governs  whether  a  minimum,  maximum,  or  specific  mean  response  is 
required.  Combining  the  qualities  of  accuracy  and  robustness  can  strengthen  the 
understanding  of  each  RPD  approach  and  a  more  informed  evaluation  of  the  competing 
procedures  can  be  made. 

The  proposed  approach  for  using  accuracy  and  robustness  to  compare  multiple 
RPD  strategies  is  now  presented  for  a  general  problem  in  which  there  exist  K  strategies 
and /responses.  For k  =  1,2,..., AT  and  j  =  1,2,...,/,  define  the  following: 

•  xk  is  the  robust  point  found  using  RPD  strategy  k 

•  m[pJ  and  v(k'’j  are  the  system’s  predicted  mean  and  variance,  respectively,  at  xk  for 
response  j 

•  m['  \ and  are  the  system’s  realized  mean  and  variance,  respectively,  atx/;  for 
response  j 

•  tj  is  the  desired  target  value  for  response  j  if  it  is  a  target-is-best  experimental 
scenario 

These  values  can  be  visualized  in  tabular  form  as  shown  in  Table  20.  The  last  row  of 
Table  20  is  labeled  “Ideal.”  The  ideal  values  are  the  best  possible  means  and  variances 
that  can  be  attained  for  each  response  of  the  system.  These  can  be  gained  analytically  or 
through  experimentation  when  known  functions  or  simulations  are  being  used.  Let  yj 

represent  the  true  response  j.  The  ideal  mean  for  response  j  is  then  defined  as: 
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[min  E[vj] 


rrij  =  t 


J 

max 


if  response  j  is  a  "smaller-the-better"  scenario 
if  response  j  is  a  "target-is-besf '  scenario 
if  response  j  is  a  "larger-the-better"  scenario 


(74) 


The  ideal  variance  for  response  j  is  defined  as 

v*  =  min  Var^yjJ  . 
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Table  20.  Tabular  Visualization  of  RPD  Results  for  K  Strategies  and  /  Responses 
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5.2.1  Accuracy  Quality  Index 

Previously,  the  accuracy  of  a  RPD  strategy  was  defined  as  how  close  the 
predictions  of  the  system  are  to  the  actual  realizations  of  the  system  when  it  is  continually 
executed  at  a  specific  operating  point.  Therefore,  the  row  in  Table  20  labeled  “Ideal”  can 
be  ignored  at  this  time.  Since  the  values  in  Table  20  are  most  likely  in  different  units  or 
scales,  each  column  of  the  table  containing  the  “Predicted”  and  “Realized”  rows  must 
first  be  normalized  to  the  range  [0,1]  to  eliminate  these  effects.  Nonnalized  values  are 
denoted  with  For  example,  m).p’  or  v[r\ .  Now  define  the  following  vectors: 
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•  pA  =  \jh[Pi  ,v[p?  is  the  vector  of  normalized  predicted  means  and 

variances  for  RPD  strategy  k 

•  rk  =  \jh(kl,v(k},...,m ([-'), v[7j]is  the  vector  of  nonnalized  realized  means  and 
variances  for  RPD  strategy  k 

The  Accuracy  Quality  Index  for  strategy  k  across  all  /  responses,  denoted  Ak ,  is 
defined  as  a  measure  of  the  distance  between  its  vector  of  nonnalized  predicted  values  pA 
and  its  vector  of  normalized  realized  values  rk .  It  is  computed  using  Equation  (75). 

4=H|P.-r<ll/^7  (76) 

The  distance  measure  in  Equation  (76)  is  first  scaled  by  the  maximum  possible  distance 
between  any  two  points  in  the  2.7-dimensional  unit  hypercube  and  is  then  subtracted  from 
1.  This  procedure  results  in  Ak  being  contained  in  the  interval  [0,1]  in  which  larger  values 
are  prefened. 

Accuracy  can  be  illustrated  using  the  single  response  synthetic  case  study  from 
Section  3.3.1.  To  recap,  the  combined-array  RSM  approach  using  the  Std,  NN,  CNN,  and 
CCN  models  were  tested  against  the  combined-array  MDM  approach  using  the  KR  and 
RBF  models.  The  predicted  and  realized  means  and  variances  for  each  of  the  6  RPD 
strategies  are  shown  in  Table  9.  Figure  17  plots  the  normalized  predicted  means  and 
variances  for  each  RPD  procedure  against  their  nonnalized  realized  values.  The  accuracy 
is  a  measure  of  the  length  of  the  dashed  lines  connecting  the  predicted  and  realized 
values.  It  is  obvious  that  the  MDM  approach  using  the  RBFNN  model  is  the  most 
accurate  whereas  the  RSM  approach  using  the  standard  model  is  the  least  accurate. 
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Figure  17.  Accuracy  of  Each  Strategy  for  the  Single  Response  Synthetic  Case 
5.2.2  Robustness  Quality  Index 

The  robustness  of  a  RPD  strategy  relative  to  its  competitors,  on  the  other  hand,  is 
defined  as  how  close  the  system’s  actual  realizations  are  to  the  desired  mean  and  variance 
targets.  Thus,  the  rows  in  Table  20  labeled  “Predicted”  can  now  be  ignored.  Again,  each 
column  of  the  table  containing  the  “Realized”  and  “Ideal”  rows  must  be  normalized  to  the 
range  [0,1].  Now,  with  once  more  denoting  normalized  values,  define  the  following 
vectors: 

•  rk  =  \Jhkl  ,vk}  ,vk] ]  is  the  vector  of  nonnalized  realized  means  and 

variances  for  RPD  strategy  k 

•  I*  =  [ra*,v*,...,m*,v[]is  the  vector  of  normalized  ideal  means  and  variances 

It  should  be  noted  that/*  may  not  represent  a  combination  of  system  means  and  variances 
that  can  actually  be  achieved  by  any  control  factor  setting  in  the  design  space.  This  is 
why/*  is  treated  as  an  ideal,  or  utopian,  vector  of  system  means  and  variances. 
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The  Robustness  Quality  Index  for  strategy  k  across  all  /  responses,  denoted  Rk ,  is 


defined  as  a  measure  of  the  distance  between  its  vector  of  nonnalized  realized  values  rk 


and  the  vector  of  normalized  ideal  values/* .  It  is  computed  using  Equation  (76). 


**=1- 


rk~I 


lj2J 


(77) 


Similar  to  Ak ,  larger  values  of  Rk  in  the  interval  [0,1]  are  preferred.  The  synthetic  case 

study  from  Section  3.3.1  can  further  illustrate  the  robustness  quality.  The  goal  of  the 
study  was  to  locate  the  operating  point  that  resulted  in  a  response  of  8  with  minimal 
variance.  In  this  case,  the  smallest  achievable  variance  in  the  design  space  is  18.  This  is 
found  by  minimizing  Equation  (54).  Therefore,  the  ideal  vector  was  I*  =  [8, 1 8] .  Figure 

18  plots  the  nonnalized  realized  means  and  variances  for  each  RPD  procedure  against 
their  nonnalized  ideal  values.  Robustness  is  a  measure  of  the  length  of  the  dashed  lines 
connecting  the  realized  and  ideal  values.  The  operating  point  found  via  the  MDM 
approach  using  the  RBFNN  model  stands  out  as  the  most  robust  of  the  6  solutions. 

As  a  matter  of  fact,  Figure  18  also  shows  that  points  labeled  “cT  are  solutions  that 
are  completely  dominated  by  another  solution.  A  solution  is  dominated  if  another 
solution  resulted  in  a  better  realized  mean  and  a  better  realized  variance.  In  this  case,  Std 
is  dominated  by  the  other  5  solutions.  Also,  NN  and  CNN  are  dominated  by  CCN,  KR, 
and  RBF.  CCN  is  dominated  by  KR  and  RBF.  Finally,  KR  and  RBF  are  considered  non- 
dominated  solutions. 
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Figure  18.  Robustness  of  Each  Strategy  for  the  Single  Response  Synthetic  Case 
5.2.3  Joint  Quality  Index 

The  perfect  RPD  strategy  has  two  characteristics.  First,  it  precisely  predicts  the 
realized  means  and  variances  for  each  response.  Second,  it  realizes  the  ideal  means  and 
variances.  These  characteristics  equate  to  having  an  Accuracy  Index  and  a  Robustness 
Index  equal  to  1 .  The  accuracy  and  robustness  indices  can  now  be  used  to  determine  the 
overall  Joint  Quality  Index  for  each  strategy  k  by  utilizing  the  equation 

a=l-||q,-l||/V2  (78) 

where  qA  =  [  Ak ,  Rk  ]  and  1  =  [  1 ,1  ] .  Again,  the  distance  measure  in  Equation  (78)  is  scaled 
by  the  maximum  possible  distance  between  any  two  points  in  the  unit  square  and  is  then 
subtracted  from  1 .  The  strategy  resulting  in  the  largest  value  for  Qk  is  then  determined  to 

be  the  best  complete  procedure — out  of  the  K  competing  strategies — for  solving  the  RPD 
problem  across  the  J  responses  concurrently.  Figure  19  plots  the  Accuracy  Index  and  the 
Robustness  Index  for  each  RPD  strategy  in  the  synthetic  case  study  from  Section  3.3.1. 
The  overall  joint  quality  index  is  simply  a  measure  of  the  distance  between  each 
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strategy’s  accuracy-robustness  pairing  and  the  point  (1,1).  Again,  the  operating  point 
found  via  the  MDM  approach  using  the  RBFNN  model  shows  to  have  the  best  overall 
quality. 
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Figure  19.  Overall  Joint  Quality  of  Each  Strategy  for  the  Single  Response  Synthetic 

Case 

The  accuracy,  robustness,  and  overall  quality  indices  that  correspond  to  Figure  17, 
Figure  18,  and  Figure  19  are  displayed  in  Table  21 .  It  is  clear  that  for  this  case  study,  of 
the  6  RPD  procedures,  the  MDM  approach  that  utilizes  the  RBF  model  is  superior.  A 
comment  must  also  be  made  in  regards  to  the  Std  model’s  robustness  score  of  0.  The 
RSM  approach  that  utilizes  the  standard  response  surface  model  has  been  proven  to  be  a 
very  successful  strategy.  However,  the  accuracy,  robustness,  and  overall  quality  indices 
as  defined  here  are  relative  to  the  K  strategies  being  compared  in  the  study.  Therefore,  as 
it  relates  to  the  other  five  competitors  in  this  particular  problem,  the  RSM  approach  using 
the  Std  model  does  not  produce  a  very  robust  solution. 
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Table  21.  Quality  Indices  for  Each  Strategy  in  the  Single  Response  Synthetic  Case 


RPD 

Strategy 

k 

Accuracy 

Ak 

Robustness 

Rk 

Overall  Quality 
Qk 

RSM  w /  Std 

1 

0.18  (6) 

0.00  ( 6) 

0.09  ( 6) 

RSM  w /  NN 

2 

0.62  ( 4 ) 

0.61  (4) 

0.62  (4) 

RSM  w /  CNN 

3 

0.62  (4) 

0.61  ( 4 ) 

0.62  (4) 

RSM  w /  CCN 

4 

0.76  (3) 

0.67  (3) 

0.71  (3) 

MDM  w/  KR 

5 

0.79  (2) 

0.72  (2) 

0.75  (2) 

MDM  w /  RBF 

6 

0.94  (7) 

0.75  (7) 

0.82  (7) 

Rank  in  Parenthesis  ( ) 


5.3  Application  and  Results 

The  benefits  of  a  more  holistic  RPD  strategy  assessment  using  the  quality  index 
will  be  demonstrated  with  a  discrete  event  simulation.  The  ND  procedure  from  Chapter 
IV  was  assessed  against  the  KY procedure  using  six  different  RPD  strategies:  the 
combined-array  RSM  approach  that  employs  the  Std,  NN,  CNN,  and  CCN  models  and  the 
MDM  approach  that  employs  the  KR  and  RBF  models.  The  proposed  quality  indices 
were  then  used  to  compare  the  12  different  multi-response  strategies  listed  in  Table  22. 


Table  22.  Multi-Response  RPD  Strategies  Used  for  Comparison 


Label 

Strategy 

Modeling  Approach 

j^yStd 

Minimize  KY 

RSM  w/  Standard  Model 

NDs,d 

Maximize  D 

RSM  w/  Standard  Model 

kynn 

Minimize  KY 

RSM  w/  Noise-by-Noise  Model 

NDnn 

Maximize  D 

RSM  w/  Noise-by-Noise  Model 

kycnn 

Minimize  KY 

RSM  w/  Control-by-Noise-by-Noise  Model 

NDcnn 

Maximize  D 

RSM  w/  Control-by-Noise-by-Noise  Model 

j^yCCN 

Minimize  KY 

RSM  w/  Control-by-Control-by-Noise  Model 

ndccn 

Maximize  D 

RSM  w/  Control-by-Control  -by-Noise  Model 

KYkr 

Minimize  KY 

MDM  w/  Kriging  Model 

ndkr 

Maximize  D 

MDM  w/  Kriging  Model 

j^yRBF 

Minimize  KY 

MDM  w /  RBFNN  Model 

ndrbf 

Maximize  D 

MDM  w /  RBFNN  Model 
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A  multi-response  RPD  study  was  performed  using  Kelton  et  al.’s  [84]  automotive 
maintenance  and  repair  shop  (AMRS)  simulation  (Model  6-1)  as  it  was  developed  in 
Arena®.  The  two  performance  measures  of  interest  were  the  average  daily  profit  and  the 
average  daily  number  of  late  wait  jobs,  labeled  y1  and y,  respectively.  The  responses  were 
averages  of  100  independent  replicates  at  a  specified  design  setting.  The  input  factors 
that  influence  Vj  andy2  are  listed  in  Table  23.  The  factors  xl  and  x2  were  controllable  within 
their  minimum  and  maximum  values.  Also,  since  the  desire  was  to  find  the  control  factor 
settings  that  were  robust  to  uncertain  demand,  the  number  of  calls  that  arrive  to  the  shop 
each  day  was  a  Poisson  random  variable  with  meanzj .  The  mean  of  the  Poisson 

distribution  itself  was  assumed  to  be  normally  distributed  with  a  known  mean  and 
standard  deviation.  This  is  comparable  to  the  approach  taken  by  Wild  and  Pignatiello 
[44]  in  which  they  found  a  job  shop  design  that  was  robust  to  uncontrollable 
environmental  factors  such  as  the  mean  inter- arrival  times  for  parts.  Finally,  in  order  to 
structure  the  RBFNN  metamodel  so  it  was  well-generalized  with  a  minimal  risk  of  over¬ 
fitting  the  experimental  data,  10  rounds  of  the  10-fold  cross  validation  procedure  were 
used. 


Table  23.  Input  Factors  for  the  AMRS  Problem 


Label 

Factor 

Min 

Max  Mean 

Std 

Dev 

Xi 

Maximum  work  hours  available  per  day 

20 

40 

x2 

Service  buffer  hours  allowed  for  waiting 
customers 

0.5 

2 

zi 

Mean  number  of  calls  per  day 

25 

3 
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The  goal  of  this  RPD  study  was  to  identify  the  robust  setting  o  f  x,  andx,  that 
jointly  maximized  the  average  daily  profit  (a  larger-the-better  case)  and  minimized  the 
average  daily  number  of  late  wait  jobs  (a  smaller-the-better  case)  with  minimum 
variability  for  each  response.  Each  response  was  treated  with  equal  weighting.  Each 

-5 

mean  and  its  associated  variance  were  also  treated  with  equal  weighting.  A  5  full- 
factorial  combined  array  design  was  used  to  build  each  model.  The  mean  and  variance 
models  for  the  Std,  NN,  CNN,  CCN,  KR,  and  RBF  response  models  can  be  viewed  in 
Appendix  B.  Table  24  summarizes  the  predicted,  realized,  and  ideal  means  and  variances 
ofy^  andy2  at  the  12  robust  points.  The  realized  values  are  the  results  of  100  simulations 
of  each  robust  operating  point  using  common  random  numbers  for  the  noise  factor. 
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Table  24.  Multi-Response  RPD  Results  for  the  AMRS  Simulation 


RPD 

Robust  Point 

Mi 

of 

Mi 

oi 

Strategy 

*2 

Goal 

Max 

Min 

Min 

Min 

0.21 

1.00 

Predicted 

576.68 

637.39 

0.43 

8.03  XI 0'4 

Realized 

565.11 

375.86 

0.43 

4.05  X10'4 

NDStd 

-0.10 

1.00 

Predicted 

561.95 

601.28 

0.41 

7.56  X10'4 

Realized 

578.70 

224.51 

0.39 

2.86  X10'4 

j£Ym 

0.21 

1.00 

Predicted 

576.60 

641.85 

0.42 

1 .35  X  HP 

Realized 

565.11 

375.86 

0.43 

4.05  X10'4 

-0.10 

1.00 

Predicted 

559.87 

605.74 

0.40 

1.30X  10'3 

XI) 

Realized 

578.70 

224.51 

0.39 

2.86  X10'4 

kycnn 

0.20 

1.00 

Predicted 

574.47 

641.76 

0.42 

1.60  xi  O'3 

Realized 

566.22 

349.73 

0.43 

5.58  X10'4 

NDcnn 

-0.18 

1.00 

Predicted 

Realized 

550.86 

573.68 

599.38 

172.56 

0.39 

0.38 

1.08  X  10'J 
3.81  X10'4 

j£Yccn 

0.20 

1.00 

Predicted 

574.47 

563.20 

0.42 

1.56  x  10"J 

Realized 

566.22 

349.73 

0.43 

5.58  X10'4 

0.14 

1.00 

Predicted 

573.85 

566.06 

0.42 

1 .50  x  or 

ND 

Realized 

570.48 

410.09 

0.42 

3.01  X10'4 

KYkr 

-0.06 

0.83 

Predicted 

583.46 

52.51 

0.43 

6.44  XI  O'3 

Realized 

581.85 

357.48 

0.44 

2.54  X10'4 

-0.29 

1.00 

Predicted 

565.37 

37.58 

0.36 

2.47  XI  O'3 

ndkr 

Realized 

561.05 

124.41 

0.37 

3.63  XI  O'4 

-0.04 

1.00 

Predicted 

570.77 

122.08 

0.39 

5. 60  XI 0'4 

Realized 

576.11 

304.20 

0.40 

2.91  X10'4 

ndrbf 

-0.12 

1.00 

Predicted 

569.78 

122.48 

0.38 

5. 47  XI 0'4 

Realized 

580.86 

209.28 

0.39 

2.43  XI  O'4 

Ideal 

582.04 

16.70 

0.28 

1.59  X10'4 

The  first  thing  to  note  while  examining  Table  24  is  that  it  is  very  cumbersome. 
This  study  compares  the  results  of  12  procedures  across  only  two  responses.  It  is  difficult 
to  inspect  this  table  and  make  strong  conclusions  regarding  the  study.  However,  utilizing 
the  accuracy,  robustness,  and  joint  quality  measures  allows  for  a  more  comprehensive 
assessment  of  the  competing  strategies.  Illustrations  of  these  measures  for  each  response 
are  shown  in  Figure  20  while  the  indices  across  both  responses  are  displayed  in  Table  25. 
Based  on  their  quality  indices,  it  is  clear  that  NDRBF  and  KYRBF  are  superior  to  the  other 
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10  strategies.  This  has  two  implications.  First,  the  mean  and  variance  predictions  at  their 
robust  points  are  comparable  to  how  the  system  actually  performs  at  their  robust  points. 
This  implies  that  the  mean  and  variance  models  generated  via  the  MDM  approach  using  a 
RBFNN  model  are  accurate.  Second,  the  system’s  actual  perfonnance  across  the  two 
responses  at  their  robust  points  is  more  robust  than  the  other  strategies.  This  illustrates 
that,  in  this  case,  utilizing  the  MDM  approach  with  a  RBFNN  model  for  an  RPD  study  of 
a  simulation  is  superior  to  the  RSM  approach  that  uses  the  polynomial  regression  models. 

It  is  also  interesting  to  note  the  comparison  between  the  ND  procedure  and  the  KY 
procedure.  The  ND  procedure  outperfonns  the  KY  procedure  for  every  modeling  strategy 
in  terms  of  robustness.  That  is,  R2  >  Rn  R4>R,,  R6>Rs,  Rs>R7 ,  Rw  >  R9 ,  and 
Rl2  >  Ru  ■  Five  of  the  top  six  approaches  in  terms  of  robustness  utilize  the  ND  procedure. 

Also,  NDkr  and  NDRBF  are  the  only  non-dominated  solutions  across  both  responses. 
Similarly,  the  ND  procedure  surpasses  the  KY  procedure  for  every  modeling  strategy  in 
terms  of  overall  joint  quality.  This  further  demonstrates  that  better  joint  robust  points  can 
be  located  by  first  putting  the  means  and  variances  of  all  responses  on  a  level  playing 
field. 
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Accuracy 


Robustness 


Overall  Joint  Quality 


Predicted  Realized 

KY3*1  □  ■ 

NDNS  o  • 

KY □  ■ 

M>w  O  • 

JCJOW  □  ■ 

JVD®*  O  • 

jn«c,v  □  ■ 

NDCCS'  O  • 

A3**  □  ■ 

JVD**  O  • 

ATJKW  □  ■ 

ndrbf  o  • 

Ideal  x 


Figure  20.  Accuracy,  Robustness,  and  Overall  Joint  Quality  of  Each  Strategy  for 

the  AMRS  Simulation 
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Table  25.  Quality  Indices  of  Each  RPD  Strategy  for  the  AMRS  Simulation 


RPD  Strategy 

k 

A 

k 

R 

k 

Q 

k 

j^yStd 

1 

0.71 

(7) 

0.17 

(9) 

0.38 

(9) 

NDStd 

2 

0.58 

(9) 

0.52 

(2) 

0.55 

(3) 

KYm 

3 

0.73 

(6) 

0.17 

(9) 

0.38 

(8) 

NDm 

4 

0.56 

(10) 

0.52 

(2) 

0.54 

(4) 

kycnn 

5 

0.71 

(8) 

0.11 

(11) 

0.34 

(12) 

NDcnn 

6 

0.49 

(11) 

0.49 

(4) 

0.49 

(6) 

J^yCCN 

7 

0.77 

(5) 

0.11 

(11) 

0.35 

(11) 

NDccn 

8 

0.83 

(1) 

0.26 

(8) 

0.46 

(7) 

KY ™ 

9 

0.44 

(12) 

0.33 

(7) 

0.38 

(10) 

ND™ 

10 

0.80 

(4) 

0.36 

(6) 

0.52 

(5) 

KYrbf 

11 

0.82 

(2) 

0.43 

(5) 

0.58 

(2) 

ndrbf 

12 

0.81 

(3) 

0.56 

(1) 

0.66 

(1) 

Rank  in  Parenthesis  ( ) 


5.4  Summary 

Chapter  V  presented  a  comprehensive  assessment  framework  for  comparing  the 
results  of  several  RPD  problem  solving  strategies.  This  approach  expanded  beyond  the 
typical  tactic  of  simply  investigating  how  close  the  system’s  predicted  means  and 
variances  at  each  of  the  established  robust  points  are  to  ideal  target  values.  When  it  is 
feasible,  the  use  of  confirmation  experiments  permits  further  examination  into  the 
accuracy  and  robustness  qualities  of  each  method.  This  approach  allows  for  a  more 
knowledgeable  evaluation  of  the  competing  procedures. 
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VI.  Conclusions 


This  dissertation  addressed  three  fundamental  objectives.  The  first  objective  was 
to  broaden  the  combined-array  RDM  approach  that  relies  exclusively  on  low-order 
polynomial  models.  The  second  objective  was  to  develop  an  approach  for  multi-response 
RPD  problems  that  provides  a  collaborative  solution  that  is  balanced  across  the  means 
and  variances  of  each  response.  Finally,  the  third  objective  was  to  generate  a  framework 
for  evaluating  competing  RPD  problem  solving  strategies.  Chapters  III-V  detailed  the 
methodology  for  achieving  these  objectives. 

6.1  Original  Contributions 

Chapter  III  extended  the  combined-array  RSM  approach  to  include  the 
application  of  non-linear  metamodels.  The  proposed  MDM  approach  replaced  the  low- 
order  polynomial  models  with  Kriging  and  RBFNN  models.  Then,  via  the  Multivariate 
Delta  Method,  mean  and  variance  models  were  generated  from  second-order  Taylor 
series  approximations  of  the  Kriging  and  RBFNN  models.  Finally,  an  existing 
optimization  problem  that  employed  these  approximations  was  solved  to  identify  the 
robust  control  parameter  setting.  The  combined-array  MDM  approach  was  compared 
with  two  current  RPD  strategies.  First,  when  compared  to  the  combined-array  RSM 
approach  that  uses  polynomial  models,  the  MDM  approach  demonstrated  improved 
predictive  models  of  a  computer  simulation’s  mean  and  variance.  Second,  the  MDM 
approach  was  shown  to  generate  results  that  are  approximately  equivalent  to  the 
stochastic  emulator  approach  at  a  significantly  reduced  computational  cost. 
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Chapter  IV  proposed  a  multi-response  RPD  procedure  based  on  well-known 
desirability  functions  that  showed  two  benefits.  First,  it  placed  each  response,  as  well  as 
their  means  and  variances,  on  common  ground.  This  increased  the  opportunity  to 
identify  a  solution  that  is  well-balanced  across  the  means  and  variances  of  each  response. 
Second,  it  allowed  a  decision  maker  to  state  their  personal  preferences  for  the  responses’ 
means  and  variances.  The  resultant  operating  point  is  a  system  setting  that,  whether  one 
is  examining  the  problem  from  a  MSE  POV  or  a  desirability  POV,  produces  a  mutually 
robust  set  of  responses  the  decision  maker  considers  acceptable. 

Chapter  V  describes  a  framework  for  comparing  and  contrasting  competing  RPD 
problem  solving  strategies  via  several  quality  measures.  This  approach  expanded  beyond 
the  typical  tactic  of  only  investigating  how  close  the  system’s  predicted  means  and 
variances  at  each  of  the  established  robust  points  are  to  ideal  target  values.  By 
perfonning  confirmation  experiments,  further  insight  can  be  gained  into  each  RPD 
approach’s  accuracy  and  robustness  qualities.  Accuracy  measures  of  how  close  the 
predicted  mean  and  variance  of  the  system  are  to  the  realized  mean  and  variance  of  the 
system  when  it  is  repeatedly  executed  at  a  specific  operating  point.  Robustness,  on  the 
other  hand,  assesses  how  close  the  system’s  realized  mean  and  variance  are  to  their 
desired  target  values.  The  combination  of  accuracy  and  robustness  can  increase  an 
analyst’s  understanding  of  each  RPD  approach  and  allow  them  to  make  a  more  informed 
evaluation  of  the  competing  procedures. 
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6.2  Future  Research 


This  research  focused  on  the  utilization  of  Kriging  and  RBFNN  metamodels. 
However,  there  are  numerous  other  non-linear  modeling  techniques.  There  is  an  obvious 
opportunity  to  extend  this  research  to  examine  the  effect  other  metamodeling  techniques, 
such  as  spline  regression  or  other  neural  networks,  have  on  the  overall  robust  solution 
quality. 

A  second  opportunity  for  research  focuses  on  the  generalization  of  radial  basis 
function  neural  networks.  If  too  many  neurons  are  used,  then  the  overall  generalization 
of  the  network  will  be  deficient.  On  the  other  hand,  if  too  few  neurons  are  used,  the 
network  will  not  be  able  to  sufficiently  leam  the  training  data.  This  research  employed  a 
cross-validation  procedure  to  determine  the  structure  of  the  RBFNN  so  that  the  resulting 
function  was  well-generalized  with  minimal  risk  of  over-fitting  the  experimental  data. 
Opportunities  exist  for  further  heuristic  development  that  identifies  a  “robust  neural 
network  structure. 

A  third  area  of  future  research  is  concerned  with  how  critics  view  ANNs  in 

general.  This  dissertation  showed  that,  by  using  RBFNNs,  a  robust  solution  could  be 

generated  that  is  as  good  as,  and  in  some  cases  significantly  better  than,  those  produced 

from  strategies  using  standard  polynomial  regression  models.  However,  critics  of  neural 

networks  typically  have  points  of  view  similar  to  that  of  Myers  and  Montgomery: 

Our  view  is  that  neural  networks  are  a  complement  to  the  familiar 
statistical  tools  of  regression  analysis,  RSM,  and  designed  experiments, 
but  certainly  not  a  replacement  for  them,  because  a  neural  network  can  at 
best  only  give  a  prediction  model  and  not  fundamental  insight  into  the 
underlying  process  mechanism  that  produced  the  data.  [8] 
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Although  there  is  already  abundant  literature  on  identifying  salient  features  using  ANNs, 
perhaps  there  is  some  potential  for  translating  the  weights  and  interconnected  neurons  of 
an  ANN  to  the  standard  regression  coefficients  that  the  experimental  community  is  most 
familiar  with. 
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Appendix  A.  Mean  and  Variance  Derivations  for  a  Function  of  a  Random  Variable 


A.l  Univariate  Case 

Let  Y  =  f  ( X )  be  a  function  of  the  normally  distributed  random  variable  X  where 
E[X ]  =  jux and Var(X)  =  crx  .  Also,  let£  ~  a(0, cL1 ) .  Finally,  X is  independent  of  e.  If 

Y  =  f{X)  =  f(jux  )  +  f'(jux  )(X  -jux)  +  \f\vx  \X  -  vx  )2  +  ^  (79) 

is  a  second-order  Taylor  series  approximation  of  Y  centered  at  the  point  a  =  /ux ,  then  an 
estimate  for  the  mean  of  Y  is 

£[r]  =  £[/(*/,)  +  fiMxXX  -  /'*)  +  {/"(/'*)(*  -  Vxf  +  *] 

=  f(h)  +  /  '</'■,)£[£  -  +  -  /<,r]  +  i'H  (80) 

Similarly,  an  estimate  for  the  variance  of  7  is 

('«-[)']  -  I  ar[  /(/(,  )+  /'(//,  )(£-/;,)-(  /•(«,  )t  X  -u, )’  +£) 

-  !>«'[  /'!/(.,  )(  \'  -  /(,  )■  (  / "l//,  )(Y-/I,  )2]  +  Var[s] 

= £[(/'(//, x :-*•  -  ft) + !/■(/<* xx  -  nx y  y 

-E[fXfix)(X-»x)+\f\fix)(X-vsyy+c ,]  (81) 

=fx^xyE[{x-^xf]+r{^x)r^x)E[{x-MxY] 

+ \f(MS  y  e  [(x -  iax  )4  ] . -  [i  f(Mx  k  ]! + 

^nflxy-cr2x+if'(f‘xy<+^ 
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The  following  central  moments  for  X  were  used  to  derive  Equations  (80)  and  (81): 


E[X-jux  ]  =  0 

(82) 

E[(X-jux)2]  =  *2x 

(83) 

E[(X-jUxf]  =  0 

(84) 

V 

cn 

II 

■'t 

* 

A 

1 

V 

(85) 

A.2  Multivariate  Case 

Let  Y  =  f  (X)  be  a  function  of  the  normally  distributed  random  vector  X  where 
£[X]  =  nx  and Var(X)  =  Zx  .  Also,  let ~  A(0, cr  ) .  Finally,  X  is  independent  of  s.  If 

Y  =  /(X)  =  /(px)  +  V/0ix)'(X  -  Jix)  +  Kx  -  |ix)'HGix)(X  -px)  +  £  (86) 

is  a  second-order  Taylor  series  approximation  of  Y  centered  at  the  vector  a  =  px  ,  then  an 
estimate  for  the  mean  of  Y  is 

E[Y]  =  E  [/Qix  )  +  V/  0ix)'(X  -  |ix)  +  i  (X  -  jix  )' HQix  )(X  -  |ix  )  +  *] 

=  /(^x)  +  V/(iixy£[X-px]  +  i£[(X-pxyH(px)(X-px)]  +  £^]  (87) 

=  /(Hx)  +  ifr(H(Hx)2x) 

Similarly,  an  estimate  for  the  variance  of  Y  is 
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Var  [7]  =  Var  [/(fix)  +  V/  (jix)'(X  -  |ix)  +  y  (X  -  fix)'H(fix)(X  -  Hx)  +  ^] 

=  Var  [V/(jix  )'(X  -  (ix)  +  |(X  -  fix)'H(^x)(X  -  jix)]  +  Var  [e] 

=  E  [(V/(|ix)'(X  -  |ix)  +  i(X  -  ^ix)'H(fix)(X  -  ^ix))2 
-^[V/(fix)(X-fix)  +  y(X-(ix)'H(fix)(X-(ix)]2  +<j2e 

=  e  [v/(nx)'(x  -  nxxx  -  »xy Vfdix) 

+V/(fix)'(X-fix)(X-[ix)'H(nx)'(X-nx) 

+  T  (x  -  fix)'H(nx)(X  -  jix)(X  -  nx)'H(nx)'(X  -  nx)] 
-[^r(H(^ix)Zx)]"+o-£2 

=  vfbix  )'e  [(x  -  fix)(x  -  nxy  ]  v/oix) 

+v/  (Hx)'£  [(X  -  (ix)(X  -  ^ix)'H(fix)'(X  -  (ix)] 

+  \E[X-  Hx)'H(fix)(X  -  |1X)(X  -  Hx)'H(fix)'(X  -  (ix)] 
-|^(H(fix)Ex)2+cT2 

(88) 

=  V/(fix)'ExV/(fix)  +  l[2^(H(fix)ExH(fix)Ex)  +  ^(H(tix)Ex)2 
_i^r(H(|lx)Sx)2  +  cr2 

=  V/(fix)'ExV/(fix)  +  l^(H(fix)ExH(fix)Ex)  +  fr2 

The  following  central  moments  for  X,  obtained  from  Mathai  and  Provost  [85]  and 
Brookes  [86],  were  used  to  derive  Equations  (87)  and  (88).  If  A  is  a  symmetric  matrix, 
then: 


E[X- |ix]  =  0 

(89) 

£[(X-^x)'A(X-^ix)]  =  tr(AEx) 

(90) 

E  [(X  -  ftx)  A(X  -  ftx)']  =  AEX 

(91) 
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^  [(x  -  Jlx  )(X  -  (ix )' A(X  -  J1X)]  =  £  [(X  -  I»x)' A'(X  -  |»x)(X  -  nx )']'  =  0  (92) 
£[(X-|ix)'A(X-Ms)(X-Mx)'A(X-|ix)]  =  2»-(ASxAEx)  +  »-(AEx)J  (93) 
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Appendix  B.  Mean  and  Variance  Models  for  the  Automotive  Maintenance  and 

Repair  Shop  Simulation 

The  dots  in  the  figures  represent  approximations  of  the  Automotive  Maintenance  and 
Repair  Shop  (AMRS)  simulation’s  true  response  mean  or  variance. 

B.l  Sid  Models 


Figure  21.  Mean  and  Variance  of  the  Sid  Model  for  the  AMRS  Simulation 
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B.2  AW  Models 
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Figure  22.  Mean  and  Variance  of  the  jY/V  Models  for  the  AMRS  Simulation 
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B.3  CNN  Models 


Figure  23.  Mean  and  Variance  of  the  C/V/V  Models  for  the  AMRS  Simulation 
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B.4  CCN  Models 


Figure  24.  Mean  and  Variance  of  the  CCN  Models  for  the  AMRS  Simulation 
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B.5  KR  Models 


Figure  25.  Mean  and  Variance  of  the  KR  Models  for  the  AMRS  Simulation 
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B.6  RBF  Models 
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Figure  26.  Mean  and  Variance  of  the  RBF  Models  for  the  AMRS  Simulation 
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the  proposed  approach  at  a  reduced  computational  cost.  Additionally,  a  multi-response  RPD  problem  solving  technique  based  on 
desirability  functions  is  presented  to  produce  a  solution  that  is  mutually  robust  across  all  responses.  Lastly,  quality  measures  are  developed 
to  provide  a  holistic  assessment  of  several  competing  RPD  strategies. 
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